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1.0  INTRODUCTION 


This  repeal  describes  the  devele^ment  of  the  Hypersonic  Low  Density  Analysis  (HYLDA) 
computer  code.  The  HYLDA  code  is  intended  to  support  the  design  of  complex  hypersonic 
shapes  in  low  density  flows  where  real  gas  effects  are  significant  The  focus  of  fliis  effort  was  to 
correctly  characterize  the  plasma  state  and  energy  distribution  in  the  flowfield  about  high- 
velocity,  high-altitude  hypersonic  vehicles.  To  do  this,  die  HYLDA  code  solves  the  three- 
dimensional  Navier-Stokes  equations  with  fiilly-coupled  finite  rate  air  chemistry. 

The  development  of  the  HYLDA  code  was  motivated  by  the  difficulty  of  current  wind 
tunnel  facilities  to  accurately  model  the  shock  layer  around  complex  vehicles  flying  at  very  high 
altitudes  and  very  high  Mach  numbers.  These  modeling  difficulties  arise  from  the  need  to 
accurately  model  nonequilibrium  chemistry  effects,  including  plasma  states.  The  HYLDA  code 
is  capable  of  accounting  for  the  effects  of  chemical  nonequilibrium  in  high  speed  flows, 
including  tracking  the  fiee-electron  concentrations  in  the  shock  layer.  AFWL  is  also  sponsoring 
research  directed  toward  upgrading  wind  tunnel  capabilities  for  hypersonic  flows  (ref.  1). 

The  goal  of  the  development  effort  was  to  provide  a  code  capable  of  predicting  trends  in  the 
hypersonic  vehicle  flowfield  properties  to  assist  in  the  identification  of  real  gas  effects  in  existing 
flight  data,  and  to  assist  in  vehicle  optimization.  The  result  of  the  development  effort  is  the 
HYLDA  code,  capable  of  predicting  heating  rates,  pressures,  forces,  and  moments  on  complex 
shapes  in  low-density  flow  at  hypersonic  Mach  numbers.  The  code  accounts  for  the  hypersonic 
flow  effects  of  nonequilibrium  air,  including  dissociation  and  ionization. 

The  development  of  the  HYLDA  code  occurred  in  two  phases.  Phase  I  consisted  of  the 
development  of  the  flow  model  and  computational  algorithm  required  for  the  calculation  of 
nonequilibrium  air  chemistry  flows  about  hypersonic  vehicles.  Results  of  the  Phase  I  effort  were 
reported  in  the  Phase  I  Interim  Technical  Report  (ref.  2).  Phase  n  consisted  of  combining  the 
Phase  I  results  into  the  HYLDA  computer  code.  The  results  of  the  Phase  n  developriunt  effort 
are  summarized  in  this  report 

Section  2.0  describes  features  of  the  HYLDA  code,  section  3.0  details  the  theory  for  the 
Navier-Stokes/finite  rate  chemistry  algorithm,  section  4.0  describes  the  thermochemistry  models, 
and  section  5.0  discusses  the  application  of  the  HYLDA  code  to  the  three  Government-supplied 
demonstration  cases.  Details  of  the  use  and  application  of  the  HYLDA  code  can  be  found  in  the 
Software  User’s  Manual  (ref.  3). 
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2.0  HYLDA  CODE  SUMMARY 


This  section  provides  details  of  the  Hypersonic  Low  Density  Analysis  (HYLDA),  including 
a  sununaiy  of  code  features  (section  2.1),  a  description  of  the  code  data  structure  (section  2.2), 
and  a  discussion  of  the  code  multiprocessing  capability  (section  2.3). 


2.1  Code  Features 

A  summary  of  the  HYLDA  code  features  is  included  in  table  1. 

Tablet.  HYLDA  Code  Features. 


•  Geaenl 

•  Venadle  (mnldple  optkn*) 

•  Uter-ftiendly  input 

•  Ntaneiie* 

•  Aibimfily-inoffoniied  lyttcm 

•  FledMe  boundaiy  cnndiiinn  qmrifintinn 

•  Imf&la.  Ganu.Sdd«l  line  rdjucetion 

•  bviicid  Bux  iplining 

•  RUty-coupled  Bniie  me  chemiftiy 

•  Fliynct 

•  FoUNtvicr4Uake< 

•  VUeoui  or  invudd 

•  1  J«ngi«r  Of  tnifeuknt 

•  UmI  tu.  •qoilibrimii  «ir 

•  Frits  fUe  «irchemuti>’ 

•  Deti|^ 

•  Modular 

•  FleuMe  dett  Hmclaie 

•  Dynamic  memofy 

•  VeOoriaed 

•  Mnltiproceftiog 

•  Pofttfaie 

Narratives  of  the  features  in  table  1  follow: 

General 

a.  Versatile  (multiple  options).  The  HYLDA  code  is  written  to  provide  the  user  with  a  variety 
of  modeling  options.  The  options  can  be  combined  to  provide  analysis  capability  for  a  wide 
variety  of  flow  problems. 

b.  User-friendly  Input.  The  multiple  options  available  to  the  user  are  supported  by  a  user- 
friendly  input  generation  processor,  included  as  part  of  the  HYLDA  code.  The  input 
generation  scheme  provides  the  user  with  default  values  for  each  input  quantity,  and 
provides  error  checking  and  retry  options  to  correct  inputs  entered  in  error.  The  input  file  is 
friUy-commented,  and  can  be  edited  to  create  new  input  files.  C^omplete  details  of  the  user 
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input  process  are  described  in  detail  in  the  Software  Users  Manual  (ref.  3). 


Phystes 

a.  Full  Navinr-Stokes.  The  HYLDA  code  is  based  on  the  full  Navier-Stokes  equatitms, 
coupled  with  a  one-temperature  air  finite  rate  chemistry  model. 

b.  Viscous  or  inviscid.  The  code  can  be  run  in  either  a  viscous  or  inviscid  mode.  The  code 
can  also  be  run  inviscid  fiv  a  number  of  steps  before  switching  to  a  viscous  calculation. 

c.  Laminar  or  turbulent  The  code  can  be  run  for  laminar  or  turbulent  flows.  The  turbulence 
model  used  here  is  the  Baldwin-Lomax  (ref.  4)  model,  and  is  intended  for  use  on  wall- 
bounded  turbulent  flows  only. 

d.  Ideal  gas  or  equilibrium  air.  In  addition  to  the  finite  rate  chemistry  gas  model  developed  fcx 
this  contract  the  HYLDA  code  is  also  capable  of  calculating  flows  for  ideal  gases  and  for 
equilibrium  air.  The  equilibrium  air  model  used  is  the  Taimehill  noodel  described  in  ref.  S. 

e.  Finite  rate  chemistry.  The  finite  rate  chemistry  model  included  in  the  HYLDA  code  is  a 
one-temperature  model  which  includes  a  modified  version  of  the  Paric-Menees  (ref.  6)  air 
chemistry  model,  the  transport  properties  nxxlel  of  Nicolet  (ref.  7),  and  the  species 
thermodynamics  model  of  McBride  (ref.  8).  Details  of  these  models  are  included  in  section 
4.0. 

Numerics 

a.  Arbitrarily  Transformed  Coordinate  System.  The  implicit  algorithm  is  written  in  terms  of 
generalized  curvilinear  coordinates,  allowing  body  fitted  meshes  to  be  used.  Body  fitted 
meshes  facilitate  geometry  definition  (even  for  complex  three-dimensional  geometries),  and 
also  facilitate  implementation  of  boundary  conditions. 

b.  Flexible  Boundary  Condition  Specification.  The  HYLDA  code  contains  a  full  range  of 
boundary  conditions  for  supersonic  flows,  including  solid  wall  no-slip  (adiabatic  or 
isothermal,  noncatalytic  or  catalytic),  solid  wall  free-slip,  symmetry,  inflow,  outflow,  and 
freestream  conditions.  These  boundary  conditions  are  treated  implicitly,  thereby 
ixuuntaining  the  stability  and  robustness  of  the  basic  algorithm. 

c.  Implicit  Gauss-Seidel  Line  Relaxation.  The  implicit  algorithm  is  based  on  the 
MacCc*rmack  implicit  algorithm  (ref.  9).  By  utilizing  the  unfactored  block  tridiagonal 
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matrix  structure  obtained  by  linearizing  the  implicit  governing  equation,  the  MacCormack 
in^Iicit  algorithm  does  not  exhibit  the  time  step  restrictions  found  in  approximate 
factorizadon  techniques.  The  block  tridiagonal  is  solved  using  Gauss-Seidel  (G-S)  line 
relaxadoo,  solving  a  tridiagonal  matrix  for  lines  along  a  prescribed  direction.  G-S  sweeps 
update  the  solution  in  the  remaining  two  directions.  Convergence  for  each  time  step  is 
obtained  in  approximately  two  to  three  iterations  (sweeps)  due  in  part  to  the  use  of  inviscid 
flux  vector  splitting  which  enhances  the  stability  of  the  solution  procedure. 

d.  Inviscid  Flux  Splitting.  The  HYLDA  algorithm  incorporates  inviscid  flux  splitting  to  more 
accurately  model  flow  physics  than  empirical  methods  using  smoothing  factors.  The  code 
uses  a  combination  of  the  MacCormack  (ref.  9)  and  Steger/Wanning  (ref.  10)  schemes  to 
maintain  stability  across  strong  shocks. 

e.  Fully-coupled  finite  rate  chemistry.  The  HYLDA  code  solves  the  finite  rate  air  chemistry 
flowfield  about  hypersonic  bodies  using  the  Navier-Stokes  equations  coupled  to  a  one- 
temperature  finite  rate  air  chemistry  model.  Fully-coupled  equation  sets  provide  the  most 
accurate  modeling  of  the  hypersonic  flowfield,  but^  are  limited  by  the  size  of  the  matrix 
equations  to  be  solved  as  the  number  of  gas  species  increases.  The  HYLDA  air  chemistry 
model  contains  9  species,  requiring  the  solution  of  block  13  by  13  matrices  for  three- 
dimensional  (3D)  flows. 

Design 

a.  Modular.  The  HYLDA  code  is  built  using  a  modular  top-down  structure  to  facilitate  the 
incorporation  of  new  technology  as  it  becomes  available. 

b.  Flexible  Data  Structure.  The  flexible  data  structure  uses  only  the  storage  needed  for  a 
particular  application.  Flexibility  is  achieved  using  a  pointer  mapping  system  which  allows 
code  enhancements  with  a  minimum  of  code  changes. 

c.  Dynamic  Memory.  The  HYLDA  code  supports  the  request  field  length  capability  on  Cray 
and  Cyber  computers,  and  allocates  only  the  field  length  needed  for  a  particular  application. 

(L  Vectorized.  Subroutines  in  HYLDA  have  been  vectorized  to  increase  code  efficiency. 

e.  Multiprocessing.  The  HYLDA  code  is  multiprocessor  compatible. 

f.  Portable.  The  HYLDA  code  is  written  in  standard  Fmiran  to  facilitate  porting  the  code  to 
the  many  minisuper-  and  super-computers  available.  The  Software  Users  Manual  (ref.  3) 
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details  features  of  the  code  that  are  machine-specific  so  that  the  user  can  elect  to  operate 
selected  machine-specific  features  if  desiiecL 

Additional  details  regarding  general  code  features  can  be  found  in  fee  Software  Usm 
Manual  (ref.  3).  Additional  details  of  the  physical  models  used  in  fee  code  are  included  in 
sectimi  4.0.  Additional  details  regarding  the  numerics  of  fee  HYLDA  implicit  algorithm  are 
included  in  section  3.0,  and  details  of  the  code  design  are  included  below. 


2J  Design  Features 

This  section  provides  details  of  fee  code  design  features  summarized  in  table  I  and 
described  above.  A  complete  description  of  fee  code  design  can  also  be  found  in  ref.  11. 

Traditionally,  engineering  codes  have  been  built  around  a  solution  algorithin.  When  fee 
algorithm  became  obsolete,  the  code  was  discarded  and  a  new  code  was  built  around  a  new 
algorithnL  With  the  advent  of  supercomputers  with  parallel  processing  and  vectorization  and 
micro-  and  macro-processing,  some  codes  have  been  built  around  a  particular  computer 
architecture  to  optimize  vectorization,  or  to  take  advantage  of  the  parallel  architecture.  Again, 
when  the  architecture  beccnnes  obsolete,  the  code  may  have  to  be  discarded.  Solution 
algwithms  and  computer  architectures  often  change  very  rapidly,  while  code  development  has 
traditionally  been  a  very  slow  process. 

The  current  effort  was  directed  at  improving  the  ability  of  the  code  development  process  to 
rapidly  and  efficiently  incorporate  new  advances  in  algorithms  and  computer  architectures 
without  requiring  major  code  modifications.  The  ability  to  rapidly  incorporate  solution 
algorithm  and  computer  architecture  advances  is  accomplished  through  a  modular  code  structure 
with  a  cmnpact,  flexible  data  structure.  The  ability  to  rapidly  exploit  advances  in  computer 
architecture  is  accomplished  by  requiring  new  coding  to  conform  to  strict  standards  of  Fortran, 
without  extensions,  to  ensure  feat  fee  base  version  is  readily  portable  to  any  computer 
architecture.  As  new  architectures  become  available,  a  version  of  the  HYLDA  code  can  be 
adapted  to  fee  new  architecture,  while  the  basic  version  remains  portable.  However,  high  levels 
of  vecttxization  can  still  be  achieved  even  while  adhering  to  the  basic  standards  of  Fortran. 

Code  Structure.  One  of  fee  goals  of  fee  HYLDA  code  development  effort  was  to  produce 
a  code  that  can  rapidly  and  efficiently  incorporate  changes  in  solution  algorithms  while  at  fee 
same  time  being  able  to  exploit  advances  in  computer  architectures.  The  code  structure  chosen 
to  satisfy  the  above  requirements  is  based  on  a  modular  code  design  with  a  compact,  flexible 
data  structure.  Details  of  the  conceptual  design  of  the  code  structure  employed  here  are  given  by 
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Hopcroft  (ref.  12). 

A  simplified  schematic  of  the  modular  code  design  is  shown  in  fig.  1.  The  controller  acts  as 
the  executive  driver  for  the  code,  and  contains  the  coding  to  create  and  read  input  files,  as  well 
as  to  allocate  and  deallocate  zones.  Through  the  input  file,  the  controller  controls  the  selection 
of  the  solution  algorithm  and  the  chemistry  model. 


Figure  1.  Top>D<nni  Code  Structure 

The  algorithm  section  contains  one  or  more  algorithms,  each  with  a  consistent  interface  to 
the  controller,  minimizing  the  changes  required  to  install  new  algtmthm  modules.  For  the 
current  development  process,  the  algorithm  section  contains  both  the  MacCormack  explicit  and 
implicit  algorithms  (ref.  9)  extended  to  include  finite  rate  chemistry  effects. 

The  chemistry  section  contains  routines  to  provide  the  flowfield  thermodynamics.  Again, 
all  routines  interface  the  controller,  thereby  making  it  simple  to  interchange  chemistry  routines. 
At  present,  the  chemistry  section  contains  three  modules:  (1)  ideal  gas,  (2)  equilibrium  air  (ref. 
5)  and  (3)  one-temperature  finite  rate  chemistry. 

The  finite  rate  chemistry  package  consists  of  reaction  paths,  species  thermodynamic 
properties,  and  species  and  mixture  transport  properties.  The  data  requirements  for  these  three 
modules  have  been  generalized  to  facilitate  the  rapid  interchange  of  models.  Flexibility  is  most 
important  in  the  area  of  finite  rate  chemistry,  as  the  reaction  models  can  vary  from  flow  to  flow, 
while  the  species  thermodynamic  and  transport  properties  methods  may  only  be  suitable  for 
certain  temperature  ranges  and  must  therefore  be  capable  of  being  replaced  or  supplemented. 
Details  of  the  models  included  in  the  HYLDA  code  are  found  in  section  4.0. 

Flexible  Data  Structure.  The  modular  code  stmeture  described  above  is  complemented  by 
a  compact,  flexible  data  structure  which  permits  code  modifications  to  be  incorporated  without 
changing  the  basic  data  structure  of  the  code.  A  three-dimensional  program  with  a  wide  range  of 
options,  such  as  the  HYLDA  code,  requires  dynamic  storage  for  effective  use  of  computer 
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menx)ry.  It  is  impractical  to  statically  dimension  arrays  with  large  numbers,  or  to  expect  the 
user  to  change  the  dimensions  for  each  problem.  Since  the  Fortran??  standard  does  not  support 
dynamic  arrays,  it  was  necessary  to  create  the  equivalent  of  dynamic  arrays  through  the  use  of 
pointers.  Real,  integer,  and  character  data  are  maintained  in  separate  arrays.  These  arrays 
contain  all  the  data  stored  by  the  code. 

Each  of  the  data  arrays  is  divided  into  sections.  The  length  of  each  section  depends  on  the 
size  of  the  problem.  The  code  automatically  determines  the  length  of  each  section,  and 
establishes  pointers  to  these  sections.  The  code  contains  five  data  arrays,  denoted  IQ,  RQR,  RQ, 
CQ,  and  CQR. 

The  IQ  array  contains  all  of  the  integer  data  required  by  the  code,  including  the  pointers  to 
the  sections  of  all  the  data  arrays.  The  IQ  array  has  seven  sections: 

•  Basic  data, 

•  Boundary  condition  data, 

•  Turbulence  model  data, 

•  Reaction  path  data, 

•  Species  data, 

•  Transpon  property  data, 

•  Convergence  tracking  data. 

In  addition  to  the  section  pointers,  the  "basic  data"  contains  information  pertaining  to  the 
computational  domain  and  the  solution  options  chosen  dtuing  the  input  creation  process. 

The  RQR  array  contains  real  reference  variables.  The  RQR  array  contains  the  same 
sections  as  the  IQ  array.  For  the  RQR  array,  the  "basic  data"  consists  of  the  reference  velocity, 
pressure,  temperature,  and  time  step  data. 

The  RQ  array  stores  the  flowfield  data  at  each  cell.  In  each  cell,  the  RQ  array  is  divided 
into  five  sections: 

•  Nonconservative  variables, 

•  Metrics, 

•  ThemKxlynamics, 
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•  Solution  vector, 

•  Additional  variables. 

Nonconservative  variables  include  grid  coordinates  and  time  step  for  each  cell,  in  addition  to  the 
velocity  components  and  internal  energy.  Mixture  thermodynamics  includes  pressure, 
temperature,  speed  of  sound,  viscosity,  thermal  conductivity,  and  specific  heat  ratio.  Additional 
variable  array  space  is  included  for  variables  not  used  in  all  models,  or  for  new  variables  being 
tested  for  inclusion  in  the  permanent  arrays. 

The  CQ  array  contains  character*  16  data,  while  the  CQR  array  contains  character*?!  data. 
The  format  of  these  arrays  is  the  same  as  for  the  IQ  array. 

The  arrays  described  above  constitute  the  entire  data  structure  for  the  code.  The  use  of 
these  arrays  facilitates  the  top-doAvn  code  structure  by  minimizing  the  number  of  arrays  that 
must  be  passed  in  subroutine  calls.  Additionally,  this  array  structure  eliminates  the  confusion 
caused  by  excessive  variable  names,  resulting  in  a  compact,  consistent  coding  scheme. 


2J  Multiprocessing  Capability 

The  HYLDA  code  can  be  run  in  either  a  single-processor  or  a  multi-processor  mode.  The 
single-processor  mode  has  been  written  in  standard  Fortran  for  ease  of  portability  between 
machines.  The  single-processor  code  has  available  a  number  of  machine-specific  statements 
which  can  be  activated  at  the  discretion  of  the  user.  Details  of  the  machine-specific  coding  are 
included  in  the  Software  Users  Manual  (ref.  3). 

The  multi-processor  capability  of  HYLDA  is  obtained  via  portable  macros  which  are 
expanded  into  machine-specific  coding  statements.  To  maintain  code  modularity  and  portability, 
code  parallelization  efforts  were  restricted  to  schemes  that  did  not  alter  the  algorithm  itself,  but 
only  altered  the  ordering  of  the  Gauss-Seidel  iterative  solution  process.  In  this  manner,  future 
changes  to  the  algorithm  will  not  require  reworking  of  the  multiprocessing  scheme. 

The  baseline  implicit  algorithm  for  the  code  contains  Gauss-Seidel  solution  sweeps  across 
the  computational  domain.  For  example,  the  j-direction  implicit  line  is  swept  in  the  i  and  k 
directions  via  standard  Fortran  DO  loops.  In  their  natural  ordering,  these  sweeps  are  not 
parallelizable  (because  of  the  use  of  the  most  recently  updated  values  at  i-1  and  k-1  in  the  G-S 
sweeps),  but  they  do  offer  a  significant  amount  of  parallelism  when  reorganized  into  a  "red- 
black"  iteration  scheme.  In  the  red-black  iteration  scheme,  each  diagonal  in  the  i-k  plane 
represents  independent  computations  that  can  be  parallelized,  although  successive  diagonals 
retain  dependence  on  previously  updated  values  at  i-1  and  k-1. 
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The  method  selected  to  parallelize  the  algorithm  creates  parallel  processes  that  remain 
viable  between  iteration  sweeps,  with  the  required  parameters  for  each  new  sweep  obtained  from 
shared  memory.  The  disadvantage  to  this  approach  is  that  some  processes  will  be  idled  while 
waiting  for  syncronization  points.  For  large  computational  problems,  though,  the  amount  of  idle 
time  is  not  significant.  Additional  details  regarding  the  parallelization  of  the  HYLDA  code  can 
be  found  in  ref.  13. 

Code  parallelization  occurred  in  two  steps;  (1)  code  restructuring  and  (2)  installing  parallel 
constructs.  Code  restructuring  required  reordering  the  HYLDA  loop  structure  in  the  implicit 
solver  to  take  advantage  of  parallelism,  and  removing  the  loop  structure  to  a  separate  subroutine. 
Installation  of  the  parallel  constructs  included  placing  shared  variables  in  shared  memory, 
identifying  (and  not  sharing)  local  variables,  spawning  slave  processes  and  providing  a  means 
for  scheduling  work  between  the  available  CPUs. 

Code  Restructuring.  In  the  original  HYLDA  code,  implicit  tridiagonal  solves  were 
performed  along  a  j-direcdon  line,  sweeping  along  the  k-direction  for  a  given  station  in  the  i- 
direction.  This  ordering  does  not  allow  for  any  parallelism  since  each  successive  tridiagonal 
solve  is  dependent  on  the  preceding  column.  To  take  advantage  of  the  parallelism  in  the 
algorithm,  these  computations  were  reordered  into  diagonal  wave  fronts.  Each  computation 
along  a  wave  front  is  then  independent  of  all  othere  in  the  same  diagonal. 

Separating  the  implicit  solver  loop  into  a  separate  subroutine  was  necessary  because  the 
parallel  model  being  used  here  is  based  on  a  subroutine  level  of  parallelism.  That  is,  separate 
copies  of  a  subroutine  are  established  on  "slave"  processors.  All  variables  needed  as  input  or 
output  must  be  included  in  either  the  calling  sequence  or  in  common.  Common  has  many 
advantages  for  the  parallel  environment  in  the  sharing  of  memory,  but  for  this  application  the 
use  of  common  was  severely  restricted. 

Two  options  were  available  for  the  operation  of  the  new  subroutine  -  it  could  either  perform 
one  tridiagonal  solve  and  return  to  ask  for  the  next  one,  or  it  could  solve  the  an  entire  wavefront 
iterating  down  the  diagonal  itself.  Tests  determined  that  the  second  option  was  more  efficient, 
so  the  new  subroutine  was  written  that  would  solve  an  entire  wavefront 

Installing  Parallel  Constructs.  The  work  to  be  parceled  out  consisted  of  individual 
tridiagonal  solves  along  each  wavefront  Slaves  were  to  be  created  once,  and  then  when  a 
wavefront  is  ready,  they  must  access  an  available  index  along  the  wave  front,  until  all  indices 
have  been  exhausted.  Once  a  wavefront  has  been  completed,  the  slaves  must  wait  until  the  next 
wavefront  is  ready,  with  the  master  process  initiating  the  next  wave  front  This  process  is 
performed  by  a  self-scheduling  loop  on  the  column  subscript  combined  with  a  barrier. 
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In  the  context  of  multitasking,  there  are  four  classes  of  variables: 

•  Input  constants, 

•  Input  variables  that  may  change  from  iteration  to  iteration, 

•  Shared  input/output  variables, 

•  Variables  which  are  local  to  the  current  iteration. 

Input  constants  may  be  passed  through  the  calling  sequence.  Input  variables  which  may 
change  during  program  execution,  as  well  as  shared  variables  must  be  placed  in  shared  common. 
Local  variables  do  not  need  to  be  referenced  outside  the  subroutine,  and  need  not  be  passed  in 
the  calling  sequence  or  in  shared  common. 

Two  different  sets  of  tools  were  investigated  for  installation  of  parallel  constructs:  the  Cray 
multitasking  library,  and  monitor  macros.  The  Cray  multitasking  library  is  a  set  of  subprograms 
implemented  for  parallel  computing  on  the  Cray.  They  have  also  been  implemented  on  the 
AUiant.  Monitor  macros  consist  of  a  set  of  macros  which  expand  out  to  code  that  implements 
parallel  control  at  a  somewhat  higher  level  than  the  Cray  multitasking  library.  Instead  of  the 
user  issuing  scheduling  commands,  the  user  provides  a  simple  scheduling  rule  as  well  as  a 
location  to  perform  the  scheduling,  and  the  monitor  actually  performs  the  scheduling.  The 
monitor  distributes  tasks  according  to  the  user-specified  rule. 

Monitors  are  built  from  four  monitor  primitives.  A  user  can  use  prefabricated  monitors 
from  the  monitor  tool  set,  or  build  new  monitors  using  the  primitives.  Porting  the  monitors  to 
other  architectures  then  requires  implementing  the  four  monitor  primitives  on  the  new 
architecture.  Since  the  monitors  are  macros,  they  require  expansion.  The  current 
implementation  of  the  monitors  is  based  on  the  macro  preprocessor,  m4.  The  m4  macro 
preprocessor  is  a  standard  Unix  utility. 

Ihe  simplest  monitor  is  GETSUB,  which  implements  self-scheduling  for  a  single  DO  loop. 
Since  the  HYLDA  code  required  more  than  a  single  DO  loop,  it  was  necessary  to  use  the 
ASKFOR  monitor,  which  offers  substantially  more  flexibility.  The  ASKFOR  monitor  allows  the 
user  to  specify  the  rule  for  assigning  tasks.  Further,  it  allows  for  control  of  a  set  of  tasks,  each  of 
which  may  be  decomposed  into  subtasks. 

HnaUy,  the  constructs  to  make  the  code  run  in  parallel  were  installed.  In  the  main  program, 
all  shared  variables  and  the  monitor  were  initialized.  In  the  implicit  solver,  it  was  necessary  to 
create  the  slaves,  initialize  the  wavefronts,  and  invoke  the  parallel  processing  subroutine.  In  the 
parallel  processing  subroutine,  it  was  necessary  only  to  execute  the  ASKFOR  monitor.  The 
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ASKFOR  monitor  then  schedules  all  work  along  the  current  wavefront 

Results.  The  parallelized  version  of  the  HYLDA  code  was  created  following  the  process 
described  above,  and  a  simple  test  case  was  calculated  to  demonstrate  the  parallelization.  The 
monitcv  macros  are  truly  portable,  with  the  same  code  running  on  Alliant,  Sequent,  and  Cray 
machines.  Results  showed  that  at  least  60%  of  the  code  was  parallelized. 

The  following  results  were  achieved  running  the  HYLDA  code  using  a  wedge  test  case  with 
a  15x5x15  grid.  Larger  grids  were  not  attempted  because  of  limited  storage  on  the  available 
computers,  and  because  of  time  constraints.  The  15x5x15  grid  provided  the  parallel  algorithm 
wavefronts  of  size  15x15,  while  solving  tridiagonals  of  length  5.  In  the  parallel  format,  this 
structure  produced  29  wave  fronts  per  Gauss-Seidel  sweep,  with  wavefronts  ranging  from  length 
1  to  15  and  then  back  to  1.  With  longer  wavefronts,  such  as  those  expected  with  more  realistic 
problems,  the  parallel  speedups  should  be  greater  than  those  demonstrated  here  for  this  small  test 
problem. 

The  data  were  collected  from  two  sources.  The  data  for  overall  performance  were  collected 
from  the  real-time  output  from  the  Unix  "time"  command,  and  are  used  to  reflect  the  CPU  time 
for  the  entire  program.  The  parallel  loop  performance  data  were  collected  from  timing 
statements  internally  inserted  before  and  after  each  update  to  the  data  set.  These  numbers  reflect 
the  wall  clock  time  required  to  execute  the  portion  of  the  code  which  ran  in  parallel.  The  ratios 
in  the  tables  below  were  computed  using  the  following  formulae: 


Speedup  = 


time  for  1  CPU 
time  for  n  CPUs 


(2-1) 


Efficiency  =  ^  100%  . 

n 


(2-2) 
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Table  Z.  Alliant  Timing  Results 


Number  of 

CPUi 

Ovenll  peffomunce 

Pmllel  loop  perfonnence 

Speedup 

Efficieacy  (%) 

Speedup 

Ef&cieocy  (%) 

1 

1.00 

100 

1.00 

100 

1.33 

77 

44 

1.92 

64 

156 

83 

Z16 

34 

3.14 

79 

Z35 

47 

3.76 

73 

Z44 

41 

67 

152 

36 

4.27 

61 

136 

32 

4.72 

59 

Tables.  Sequent  Timing  Results 


Number  of 

CPU* 

Oveiall  perfomunce 

Panllel  loop  perfonnance 

Speedup 

Efficiency  (%) 

Speedup 

Efficiency  (%) 

01 

1.00 

100 

1.00 

100 

1.60 

80 

1.86 

93 

1.99 

66 

159 

86 

126 

37 

322 

80 

147 

49 

3.80 

76 

HQH 

161 

44 

4.15 

69 

Note  that  in  table  2  above,  the  two-CPU  speedup  for  the  parallel  loop  is  less  than  1,  indicating  a 
loss  in  performance,  but  that  the  overall  performance  does  not  differ  significantly.  The  overall 
performance  reflects  the  CPU  time  spent  by  the  main  processor,  and  table  2  shows  that  the 
addition  of  the  second  CPU  did  reduce  the  amount  of  time  that  the  main  processor  spent  on  the 
problem.  The  poor  2  CPU  parallel  loop  performance  on  the  Alliant  was  caused  by  a  problem 
with  the  Alliant  operating  system  which  permitted  only  one  of  the  two  CPUs  to  be  active  at  a 
time.  Further  investigation  of  this  anomaly  was  beyond  the  scope  of  the  present  woiic.  Note, 
though,  that  the  parallel  performance  returned  to  acceptable  levels  for  the  three-CPU  case. 
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3.0  COMPUTATIONAL  ALGORITHM 


This  section  describes  the  computational  algorithm  used  by  the  HYLDA  code.  The 
algorithm  solves  the  full  Navier-Stokes  equations  using  the  MacCormack  implicit  solution 
algorithm  (ref.  9),  extended  under  this  contract  to  include  a  fully-coupled  one-temperature  finite 
rate  chemistry  nxxiel.  Solution  of  the  full  Navi^-Stokes  equations,  rather  than  reduced  forms  of 
the  Navier-Stokes  equations,  is  required  for  this  effort  since  the  primary  interest  is  in  the  physics 
of  the  viscous  shock  layer  surrounding  hypersonic  vehicles. 

The  decision  to  model  the  finite  rate  cheixtistry  effects  using  a  one-temperature  model  was 
based  on  the  fact  that  for  the  flows  to  be  considered  here  (up  to  Mach  25),  one-temperature 
noodels  have  been  shown  to  be  sufficient  to  model  the  plasma  state  around  hypersonic  vehicles, 
as  discussed  in  section  4.0.  The  increased  accuracy  afforded  by  multi-temperature  models  was 
not  significant  enough  to  offset  the  increased  computational  costs  associated  with  the  additional 
solution  equations,  therefore  a  one-temperature  coodel  was  selected  for  the  HYLDA  codb. 

The  remainder  of  this  section  details  the  HYLDA  computational  algorithm.  Section  3.1 
presents  the  governing  equations,  section  3.2  discusses  the  computational  algorithm,  and  section 
3.3  details  the  boundary  conditions  included  in  the  code. 

3.1  Governing  Equations 

Using  a  one-temperature  finite  rate  chemistry  model  allows  the  algorithm  development  to 
follow  a  general  procedure  using  the  reacting  gas  equations,  since  the  ideal  gas  and  equilibrium 
air  assumptions  are  special  cases  of  the  one-temperature  reacting  gas  case.  The  following 
algorithm  development  is  presented  in  terms  of  the  finite  rate  chemistry  algorithm.  Ideal  gas  and 
equilibrium  air  reductions  will  be  discussed  where  necessary. 

The  three-dimensional  (3D)  fiilly-coupled  Navier-Stokes/finite  rate  chemistry  equations 
consist  of  conservation  equations  for  species  mass,  x-,  y-,  and  z-momentum,  and  total  energy. 
These  equations  can  be  written  in  strong  conservation  form  in  terms  of  physical  (x,y,z) 
coordinates  as  follows: 


au'  ^  ap'  ^  ac'  ^  an' 


(3-1) 


where  U'  =  |^Pi....,p„...,pN.,pu,pv,pw,e 
variables,  =  [ coi,  ...,cu,. ,0,0, 0,0 


represents  the  vector  of  conserved  solution 
represents  the  vector  of  chemistry  source  terms. 
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and  F' ,  G',  and  H'  represent  the  x-,  y-,  and  z-direction  flux  vectors.  Here,  p*  denotes  the  density 
of  species  s,  (u,v,w)  represent  the  cartesian  components  of  velocity,  p  represents  the  total  density 
(obtained  by  summing  the  species  densities),  e  represents  the  total  energy  per  unit  volume,  and 
CO,  represents  the  species  s  production/depletion  rate,  obtained  from  the  law  of  mass  action.  In 
the  above  equations,  primes  denote  physical  (nontransformed)  vectors,  and  T  represents  a  vector 
transpose,  and  N,  denotes  the  number  of  species  present 
The  flux  vectors  F,  G',  and  H'  are  written  as: 


Pl(u-t-ttl) 


p,(u+i) 


Pn(u  +  I^) 
pu^+im„ 
pu+t,, 
puw+t„ 


3T  ^  A 

Lu(e+p>+-UT„+VT,y+WX,a-X-^+£p,U,h, 


G'= 


Pi(v  +  v,) 


Pt(v  +  v,) 


(3-2a) 


(3-2b) 


Pn.(v  +  vj^) 

PVU+^y, 

pv^+ptT^ 

PVW+Xy, 

3T  ^ 

.  v(e+p)+ux^+vT„+wx„-X,-^+2;p,v,h,. 
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Pi(w  +  w,) 


P.(w  +  W,) 


H- 


(3-2c) 


Pn(w  +  w^) 
pwu+^ 
pvw+^ 
pw*+ptT„ 

3T  ^  ^ 

.w(W-p>fUXa+Vt^+Wt„-X'^+2p,W,h,. 

oz 


N. 

where  p,  is  the  species  s  mass  density,  p=XP<  ~  density,  N,  =  total  number  of  species, 

u,v,w  »  cartesian  components  of  velocity,  w,,  »  cartesian  components  of  species  diffusion 
velocity,  p  »  pressure,  T  -  temperature,  e  >=  total  energy,  X  » thermal  conductivity,  h,  =:  species 
enthalpy,  and  tjj  is  the  shear  stress  tensor. 

For  the  ideal  gas  or  equilibrium  air  assumptions,  the  N,  species  continuity  equations  reduce 
to  a  single  equation  for  the  total  density,  and  the  di^sion  velocity  terms  (u„v„w,)  are  zero. 

The  shear  stress  matrix  is  given  in  tensor  notation  as: 


Xij=-P 


9qi  ^  3qj 

9xj  dxj 


3qi  j. 


(3-3) 


where  p  is  viscosity,  q  is  the  vector  of  canesian  velocity  components,  Pb  is  the  mixture  bulk 
viscosity  (taken  here  as  -  -|-u),  and  5  is  the  Kronekar  delta. 

The  species  diffusion  velocities  are  determined  from  (ref.  14): 


p,q,=-pD, 


ax, 

axj  ’ 


(3-4) 
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where  i  represents  any  of  the  three  diffusion  velocities,  represents  the  mixture  diffusion 
coefficient  of  species  s,  X,  denotes  the  species  molar  concentration,  and  Xj  represents  the 
coordinate  direction  corresponding  to  the  diffusion  velocity  being  calculated.  As  noted  above, 
for  ideal  gases  and  equilibrium  air,  the  diffusion  velocities  are  zero. 

For  finite  rate  air  chemistry  calculations,  the  viscosity  ()j.)  and  thermal  conductivity  (X)  are 
determined  using  transport  property  curve  fits.  Details  of  these  models  are  included  in  section 
4.0.  For  non-reacting  gases,  viscosity  is  obtained  from  Sutherland’s  law,  and  thermal 
conductivity  is  related  to  the  viscosity  using  the  Prandtl  number. 

Pressure  is  related  to  the  density  and  temperature  using  Dalton’s  Law: 

P=Zl^RT,  (3-5) 


where  N,  is  the  total  number  of  species  (for  an  ideal  gas  N,  =  1),  M,  is  the  molecular  weight  of 

A 

species  s,  and  R  is  the  universal  gas  constant 

Temperature  is  obtained  from  the  internal  energy  per  unit  mass  (ei).  The  internal  energy  per 
unit  mass  for  each  species  may  be  written  in  terms  of  its  translational  (Ci,,^).  rotational  (Ci^^), 
vibrational  (Ci^jj,),  electronic  (Ci^j^),  and  heat  of  formation  (ef )  contributions: 


Ci 


total 


~  Ci__.  +  Ci _ _  +  Cj 


‘trans 


'•rot  ^  ''•vib 


•elec 


+  er 


(3-6) 


Solving  equation  (3-6)  for  the  characteristic  temperatures  associated  with  each  form  of 
internal  energy  described  above  would  require  2*N,  +  1  equations  (assuming  the  translational 
and  rotational  temperatures  to  be  equal).  Wong  (ref.  IS)  notes  that  "...up  to  15000K,  multiple 
temperature  models  do  not  appear  to  have  demonstrated  sufficient  improvements  over  one- 
temperature  models  to  justify  the  additional  complexity  and  computational  cost"  of  using  such 
methods.  For  the  present  code  development  effort,  a  one-temperature  model  has  been  assumed. 
That  is,  all  internal  energy  modes  are  assumed  to  be  a  function  of  a  single  temperature.  Under 
the  one-tenqierature  assumption,  equation  (3-6)  can  be  rewritten  in  the  following  form: 


ei,  =  f,(T)T  +  h?, 


(3-7) 
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where  f,(T)  for  each  species  is  developed  from  curve-fit  fonnulations  for  Cp  using  a  single 
temperature,  and  h°  is  the  heat  of  formation  of  species  s.  ¥ot  ideal  gases,  the  heat  of  formation 
is  zero,  and  the  f,(T)  terms  becomes  the  specific  heat  at  constant  volume  (Cy). 

The  total  internal  energy  is  obtained  fix)m: 

N.  p,  N,  p  N.  p 

Ci  =  2  — ei.  =  2  — f.(DT  +  £  ~h?  =  fCDT  +  h?.  (3-8) 

•-1  P  ‘  .xl  P  P 


Temperature  is  obtained  from  equation  (3-8)  directly  for  ideal  gases  (with  f(T)  »  Cy),  and 
iteratively  for  a  reacting  gas.  The  equilibrium  air  model  calculates  temperature  as  a  curve-fit 
function  of  internal  energy  and  density. 

The  governing  equations  in  the  HYLDA  code  are  those  described  in  this  section  in  terms  of 
a  one-tempotiture  reacting  gas.  Equatimis  appropriate  for  ideal  gases  and  equilibrium  air  are 
subsets  of  this  general  equation  set  When  applying  the  HYLDA  code  to  an  ideal  gas  case  or  to 
a  case  using  the  equilibrium  air  irxxlel,  the  apprcqniate  governing  equations  are  obtained  thixnigh 
the  use  of  cation  flags  in  the  input  file. 

Traiuformed  Governing  Equations.  To  provide  accurate  resolution  of  flow  gradients  for 
cennplex  flowfields,  computational  grids  must  be  clustered  near  large  flow  gradients  such  as 
those  found  in  shocks,  shear  layers,  ot  in  viscous  regions  near  bodies.  Grid  clustering  results  in 
nonuniform  grid  spacings  that  are  difficult  to  incorporate  into  a  solution  algorithm  based  on  the 
physical  coordinate  system.  The  standard  approach  is  to  perform  a  generalized  coordinate 
transformation  on  the  governing  equations  to  transform  the  governing  equations  from  the 
physical  (x,y,z)  system  into  a  body-fitted  (^,11,0  system.  The  body-fitted  coordinate  system  has 
uniform  spacings  in  the  computational  domain,  facilitating  the  development  of  the  solution 
algtxithm  without  need  to  account  for  nonuniform  grid  spacings. 

The  transformed  governing  equations  are  still  written  in  strong  conservation  form; 


dt 


ac  an 
an  *  a; 


(3-9) 


whoe  the  transformed  (non-primed)  fluxes  are  defined  in  terms  of  the  physical  (primed)  fluxes 
as; 
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u=r*u'  . 


F=rH 


dx  dy  dz 


G=r* 


+  G'  +  H' 


dx 


dy 


dz 


(3-10) 


h=H 


rf 

dx  dy  dz 


w=r‘w'  , 


where  J  is  the  Jacobian  of  the  transformation,  with  primed  vectors  as  defined  in  eq.  (3-1)  and  (3- 

2). 

The  implicit  algorithm  used  by  the  HYLDA  code  is  based  on  the  governing  equations  as 
described  above.  The  next  two  sections  describe  the  develc^ment  of  the  solution  algorithm. 

3J  Algorithm  Dcvelopmait 

The  algOTithm  used  to  solve  the  governing  equations  (3-9,  3-10)  is  an  extension  of  the 
MacCormack  block  tridiagonal  implicit  algorithm  with  Gauss-Seidel  line  relaxation  and  flux 
splitting  (ref.  9).  This  contract  has  extended  the  algorithm  to  include  a  fully-coupled  one- 
temperature  finite  rate  reacting  air  chemistry  nKxlel  containing  an  arbitrary  number  of  chemical 
species. 

Solving  the  equation  set  in  fully-coupled  implicit  form  represents  die  most  physically 
orrect  representation  of  the  flow  physics,  and  should  alleviate  the  stiffness  problem  often 
associated  with  alternate  methods  of  solving  the  governing  equations  in  tlw  presence  of  finite 
rate  chemistry  (ref.  16).  The  fully-coupled  approach  will  be  computationally  intensive  if  a  large 
number  of  chemical  species  are  present.  Fot  most  applications,  the  number  of  species 
considered  is  small  Qess  than  15  species).  Considering  the  computing  power  of  current 
supercomputers,  the  increased  computing  time  required  to  solve  the  fully-coupled  system  should 
be  offset  by  the  increased  convergence  rates  obtained  by  the  fully-coupled  implicit  treatment  of 
the  equatitm  set 

The  use  of  a  fully-coupled  one-temperature  reacting  air  chemistry  model  in  HYLDA  allows 
the  solution  algorithm  to  be  written  in  general  terms  applicable  not  only  to  reacting  gas  flows. 
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but  also  to  ideal  gas  and  equilibiium  flows  without  modifying  the  basic  code.  Selection  of  the 
exact  fonn  of  the  algorithm  (ideal  gas,  equilibrium  air,  or  reacting  air)  is  accomplished  by  the 
user  through  the  input  file. 

The  solution  algorithm  will  be  developed  here  using  the  full  equation  set,  with  reductions 
for  non-reacting  gas  cases  noted  where  appropriate. 

The  implicit  algorithm  development  begins  by  recasting  the  governing  equations  (3-9),  (3- 
10)  into  delta  form.  The  first  step  is  to  split  the  momentum  fluxes  into  their  constituent  inviscid 
(Fi,Gi,Hi)  and  viscous  (Fv,Gv,Hv)  parts.  Substituting  the  split  fluxes  into  (3-9)  yields: 


au  3Fi  BGj  3Hi 

at  94  3n  ac 


aFv  aGv  aHv' 

34 

4 


(3-11) 


The  delta  form  of  the  governing  equations  if  obtained  fixMn  cq.  (3-11)  by  differentiating  with 
respect  to  t,  taking  Jacobians  with  respect  to  U  of  the  inviscid  fluxes  (Fi,Gi,Hi)  and  the  source 
term  (W),  and  by  writing  the  viscous  terms  in  terms  of  U.  Multiplying  the  result  by  At  and 
expanding  the  time  derivative  (311/30  in  a  forward  difference  yields  the  definitions: 


8ir‘  = 


IH-l 


f 


AU"  = 


(3-12) 


where  n  is  the  time  level.  In  these  definitions,  is  the  implicit  delta,  and  AU"  is  the  explicit 
delta. 

Applying  the  definitions  in  cq.  (3-12)  to  the  operations  described  above,  the  delta  form  of 
the  governing  equations  is  obtained: 
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I  +  At 


a 

Fv-4- 

n+4 

G,-4r 

N 

. 

dn 

3C 

(3-13) 


=  AU" 


where  I  is  the  identity  matrix,  A  =  9Fi/3U,  B  =  aCi/aU,  C  =  3Hi/3U  are  the  Jacobians  of  the 
inviscid  momentum  fluxes,  D  =  3W/9U  is  the  source  term  Jacobian,  and  N  is  the  matrix  relating 
the  nonconservative  vector  5V=  [5pi,..-,8ps,...,5pN,.5u,5v,5w,6T to  the  conservative 
vector  8U,  as  8V  =  (3V/3U  )8U  =  N8U.  Here,  At  is  the  time  step,  and  superscript  n  denotes  the 
time  level. 

Equation  (3-13)  represents  the  delta  form  of  the  governing  equations  for  the  solution  of  the 
Navier-Stokes  equations  with  finite  rate  chemistry. 

Flux  Splitting.  Solving  eq.  (3-13)  in  its  present  form  using  finite  difference  or  finite 
volume  techniques  would  require  some  form  of  numerical  dissipation  to  keep  the  solution  stable. 
Discrete  approximations  to  the  Navier-Stokes  equations  suffer  from  discretization  errors  which 
restrict  the  stability  and  robusmess  of  the  formulation  unless  some  form  of  artificial  dissipation  is 
added.  Numerical  dissipation  can  be  provided  by  "artificial  viscosity"  formulations,  such  as  the 
fourth  order  dissipation  function  described  by  Beam  and  Warming  (ref.  17),  or  by  a  combination 
of  second-  and  fourth-order  dissipation,  such  as  that  proposed  by  Jameson  (ref.  18).  The 
disadvantage  of  such  artificial  dissipation  functions  is  that  they  are  relatively  empirical,  and  the 
constants  can  be  adjusted  without  regard  to  the  physics  of  the  flowfield.  User-tunable  dissipation 
formulations  are  not  desirable  for  general  purpose  codes  such  as  HYLDA. 

Total  variation  diminishing  (TVD)  schemes  represent  another  method  for  adding  artificial 
dissipation  to  the  governing  equations.  TVD  schemes  have  been  shown  by  Yee  (ref.  19)  to  be 
nonoscillatory  near  strong  shocks,  but  are  more  computationally  intensive,  and  are  not  included 
in  the  HYLDA  code. 

An  alternative  form  of  dissipation  can  be  obtained  by  using  inviscid  flux  splitting,  as  first 
proposed  by  Steger  and  Warming  (ref.  10).  Initially  developed  for  the  Euler  equations,  flux 
splitting  treats  the  inviscid  terms  as  combinations  of  forward  and  backward  facing  contributions, 
depending  upon  the  sign  of  the  eigenvectors.  In  this  manner,  flux  splitting  closely  models  the 
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physics  of  the  flowfield.  Additionally,  since  the  difference  schenses  are  one-sided,  the 
formulation  is  inherently  dissipadve  (for  first  order  flux  splitting),  thereby  eliminating  the  need 
far  an  arbitrary  artificial  dissipation  term.  Flux  splitting  the  inviscid  fluxes  also  increases 
diagonal  dominance  of  the  block  tridiagonal  matrix  obtained  firom  the  implicit  governing 
equation,  thereby  enhancing  the  stability  of  the  tridiagonal  matrix  inversion. 

Since  inviscid  flux  splitting  requires  little  user-interface  to  the  model,  it  is  die  method  of 
choice  for  the  HYLDA  code.  The  flux  splitting  used  in  HYLDA  is  based  on  the  original  Steger 
and  Warming  (ref.  10)  splitting,  as  modified  by  MacCormack  (ref.  20).  A  brief  overview  of  the 
inviscid  flux  splitting  operation  is  included  here. 

Flux  splitting  is  added  to  the  implicit  (left)  side  of  the  governing  equations  (3-13)  by 
splitting  the  inviscid  flux  Jacobians  A,  B,  C  into  contributions  because  of  negative  and  positive 
eigenvalues,  as  follows: 


A  =  A'"+A~  . 


(3-14) 


where  superscript  +  denotes  positive  eigenvalue  terms  and  superscript  -  denotes  negative 
eigenvalue  terms.  The  inviscid  flux  Jacobians  B  and  C  are  also  split  as  in  eq.  (3-14). 

Flux  splitting  is  implemented  in  the  explicit  (right)  side  of  the  governing  equations  (3-13)  in 
a  slightly  different  form.  The  explicit  solution  delta  AlP  includes  inviscid  terms  of  the  form 
dFi/9^.  Steger  and  Warming  (ref.  10)  have  shown  for  an  ideal  gas  that  the  inviscid  terms  are 
homogeneous  and  of  order  one,  so  that  Fj  may  be  written  as  Fi  =  (dFi/d^  )U.  The  use  of  a  one- 
temperature  gas  chemistry  model  also  allows  the  inviscid  flux  vectors  to  retain  the  homogeneity 
property  required  for  flux  splitting,  without  requiring  the  solution  of  additional  energy  equations 
(ref.  21).  Since  the  inviscid  fluxes  are  still  homogeneous  of  order  one,  the  inviscid  flux  terms  on 
the  explicit  side  of  eq.  (3-13)  may  be  expanded  as: 


aFi 

'W 


=  ^(AU). 


(3-15) 


Equation  (3-15)  is  in  a  form  similar  to  that  described  above  for  the  implicit  side  of  the  governing 
equations,  and  again  the  flux  Jacobian  may  be  split  into  its  positive  and  negative  contributions, 
yielding: 
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(3-16) 


9Fi 


4-(Vu+a-u). 

aq 


The  operations  in  eq.  (3-15)  and  (3-16)  are  also  repeated  for  the  remaining  inviscid  fluxes  (B,  C). 
Substituting  the  results  of  eqs.  (3-15)-(3-16)  into  the  governing  equations  (3-13)  yields  the  flux 
split  form  of  the  governing  equations: 


I  +  At  + -^(T  + 

3^  84  an  an  dC  a^ 


"4 

.V.J  k  J  ^ 

-4 


8  8  ^  aHy  •  1 


(3-17) 


where  inviscid  fluxes  with  positive  eigenvalues  will  be  backward  differenced,  inviscid  fluxes 
with  negative  eigenvalues  will  be  forward  differenced,  and  viscous  terms  are  central  differenced. 
Note  that  in  eq.  (3-17),  only  diagonal  viscous  terms  are  included  on  the  implicit  Oeft)  side  of  the 
equation,  while  the  explicit  (right)  side  contains  full  explicit  terms.  Since  the  explicit  side  drives 
the  solution,  and  the  implicit  side  acts  only  as  a  convergence  accelerator,  the  final  governing 
equations  exhibit  the  character  of  the  full  Navier-Stokes  equations  with  fully-coupled  finite  rate 
chemistry. 

The  HYLDA  code  uses  two  types  of  flux  splitting  -  the  MacCormack  weak  form  and  the 
Steger/Warming  strong  form.  The  two  formulations  differ  only  in  the  indexing  selected  to 
discretize  the  inviscid  fluxes.  The  MacCormack  formulation  is  less  dissipative  than  the 
Steger/Warming  form,  and  is  applied  to  all  flow  regions  away  from  strong  shocks  and  other  large 
pressure  gradient  regions.  For  large  pressure  gradient  regions,  the  MacCormack  form  does  not 
provide  enough  dissipation  to  keep  the  solution  stable,  and  the  stronger  Steger/Warming  form  is 
applied  instead.  Switching  between  the  two  methods  is  triggered  by  the  pressure  gradient. 
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The  indexing  used  by  each  method  is  illustrated  in  figure  2  and  table  4  below.  Figure  2 
shows  a  stencil  of  two  computational  cells  in  the  i-direction.  Positive  inviscid  fluxes  travel  from 
left  to  right  in  the  figure,  while  negative  fluxes  travel  right  to  left  The  discretized  form  of  the 
governing  equations  solves  for  the  flux  across  the  cell  face  (cell  i  + 1/2). 


PosHiva  RuxaaNagativa  Fiuxas 


U1/2 

Figure!.  Flux  Splitting  StencO 


Table  4  shows  the  indexing  schemes  for  both  the  MacCormack  and  the  Steger/Warming 
formulations: 


Table  4.  Flux  Splitting  Indexing 


Positive  fluxes 

Negative  fluxes 

Method 

Implicit 

Explicit 

hnplicit 

Explicit 

MacConnack 

Steger/Warming 

ISB 

AtUi 

Aii-i 

AJiiUi,.! 

Similar  operations  apply  to  the  inviscid  fluxes  in  the  other  two  directions. 

These  flux  splitting  indices  are  used  in  the  discretized  governing  equations,  to  be  discussed 

next 

Solution  Procedure.  The  solution  procedure  for  the  governing  equations  (3-17)  consists  of 
discretization  of  the  governing  equations,  and  solution  of  the  discretized  equations.  The  flux 
split  governing  equation  (3-17)  is  discretized  using  a  finite  volume  approach,  with  negative 
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eigenvalue  inviscid  fluxes  being  forward  differenced,  positive  eigenvalue  inviscid  fluxes 
backward  differenced,  and  viscous  terms  central  differenced.  Indices  for  the  inviscid  terms  vary 
according  to  whether  MacCormack  weak  flux  splitting  or  Steger/Warming  strong  flux  splitting  is 
employed,  subject  to  the  mles  discussed  in  the  previous  section. 

The  discretized  form  of  equation  (3-17)  represents  the  fuUy-implicit  governing  equation, 
and  as  such  is  a  nonlinear  equation  for  which  practical  solution  methods  are  not  available.  To 
make  the  governing  equations  tenable,  two  operations  must  be  performed.  First,  the  nonlinear 
terms  are  linearized  to  the  previous  time  levels,  resulting  in  a  banded  block  septadiagonal 
matrix,  the  solution  of  which  is  still  not  currently  practical.  Next,  the  difference  equation  is  cast 
into  block  tridiagonal  form,  for  which  current  solution  techniques  are  practical,  by  selecting  the 
j-direction  as  the  implicit  direction.  The  deltas  in  the  i-  and  k-directions  are  relaxed  to  the 
previous  time  level  and  taken  to  the  right-hand  side,  resulting  in  a  block  tridiagonal  algorithm  of 
the  form: 

ar'SUIDii.i  +  ar'SUSl  +  af = 

AUi;j.k-a5SUl,,j.k-a45UE.i.j.k-aS5U]|j,k«+aS5USj.k.,  ,  (3-18) 


where  like  terms  have  been  collected  into  the  following  definitions: 
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ar**At-  +  . 

af*  s  At-  -Bjj-n^k  +  [  Gt,Nj  •  . 
•  ^ 

*  ^ 

.5= At-  . 

b  * 

-  ■' 

b  « 

.3= At-  +  (  Hjn]  ■  . 

. 

f  * 

•?=A.'^t,.-*(H;N]  •  . 

b  « 


(3-19) 


Here,  the  subscript  m  depends  upon  the  type  of  inviscid  flux  splitting  employed,  as  defined  in  the 
previous  section,  with  m  =  1  for  Steger/Waiming  strong  flux  splitting,  and  m  =  ‘A  for 
MacCormack  weak  flux  splitting. 

Equation  (3-18),  with  the  definitions  in  eq.  (3-19),  represent  the  discretized  governing 
equations  for  the  fully-coupled  Navier-Stokes/finite  rate  chemistry  equations. 

The  solution  technique  selected  for  the  HYLDA  code  is  based  on  the  MacCormack  implicit 
block  tridiagonal  method  using  Gauss-Seidel  line  relaxation  in  the  non-implicit  directions  (ref. 
9).  The  solution  is  updated  along  implicit  lines  using  standard  block  tridiagonal  inversion 
methods.  The  remaining  two  directions  are  updated  by  successive  solution  sweeps  through  the 
computational  domain.  In  principle,  the  iteration  sweeps  should  be  continued  until  the  deltas  in 
the  sweep  directions  converge.  In  practice,  2-3  sweeps  per  iteration  are  sufficient  unless  time 
accurate  solutions  are  required. 

The  solution  process  follows  a  one-step  implicit  update,  as  follows: 


a.  Calculate  the  explicit  delta: 


AUi;j.k=-At^ 


-§-F+-^G+^H 
A^  ATI  AC 


(3-20) 


i.j.k 


where  the  operators  and  represent  the  discretization  and  flux  splitting  operations 

described  in  the  previous  sections. 

b.  Calculate  the  implicit  delta  for  6Ujj,ic: 


i+...Uujji=Au;j.k  : 


(3-21) 


c.  Update  the  solution: 


•  (3-22) 

The  block  tridiagonal  solver  used  here  is  capable  of  solving  an  arbitrary  size  block  tridiagonal 
system  of  equations.  Thus,  no  change  in  the  solver  is  required  to  change  from  an  ideal  gas  to  a 
multiple  species  reacting  gas  case,  in  keeping  with  the  desire  to  provide  a  flexible  code  structure. 

3J  Boundary  Conditions 

The  HYLDA  code  boundary  condition  treatments  are  implemented  implicitly  to  maintain 
code  stability.  The  following  boundary  conditions  are  available: 

•  Supersonic  inflow; 

•  Supersottic  outflow; 

•  Supersonic  freestream; 

•  Solid  wall  (no-slip); 

•  Adiabatic  or  specified  temperature; 

•  Noncatalytic  or  catalytic; 
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•  Laminar  or  turbulent; 

•  Solid  wall  (free-slip); 

•  Symmetry. 

The  implicit  boundary  conditions  are  derived  from  known  explicit  boundary  conditions  by  the 
foUowing  procedure: 

a.  Specify  the  known  explicit  boundaiy  condition  in  terms  of  the  nonconservative  solution 
vector  V: 


Vbc  =  Z'V^j  +  r  .  (3-23) 

where  V=  ^  Pi,,..,p„...,pN^,u,v,w,Tj  ,  subscript  be  denotes  the  boundary  cell,  and 
subscript  adj  denotes  the  cell  across  the  boundary  surface  (adjacent)  to  the  boundary  cell. 


b.  Develop  the  implicit  boundary  condition  from  (a)  by  deriving  the  delta  form  of  the 
boundary  condition,  as  follows: 

5Vbc  =  Z'SV^j  +  V«,j5Z'  +  5r  .  (3-24) 


Now,  5V  may  be  written  in  terms  of  the  implicit  solution  vector  delta  5U  as: 


au 


5U  =  N5U  . 


(3-25) 


Substituting  eq.  (3-25)  into  eq.  (3-24)  and  solving  for  5U  yields: 


(3-26) 


8Ubc  =  [  Nte‘Z'N«,j]  5U«ij  +  [  NSl V^jSZ'  N^JSr] 
=  ZSUadj  +  X  • 


For  the  boundary  conditions  considered  here,  X  ~  ^  implicit  boundary  conditions 

currently  contained  in  the  HYLDA  code  may  be  represented  as: 


5Ubc  =  Z5U«ij 


(3-27) 


The  remainder  of  this  section  describes  the  explicit  boundary  conditions  used  in  the  code.  The 
implicit  boundary  conditions  have  all  been  derived  according  to  the  development  of  eq.  (3-27) 
above. 

Supersonic  Inflow,  Freestream.  The  supersonic  inflow  and  freestream  conditions  maintain  the 
user- specified  freestream  conditions,  as: 


Vbc  =  V..  . 


(3-28) 


where  subscript  ^  denotes  user-specified  fireestream  values. 

Supersonic  Outflow.  The  supersonic  outflow  boundary  condition  extrapolates  the  flow  in  the 
adjacent  cell  to  the  boundary  cell,  as: 


Vbc  =  V«,j  .  (3-29) 

Solid  Wall  (No-slip).  Although  the  exact  form  of  the  solid  wall  (no-slip)  boundary  condition 
depends  on  the  options  selected  by  the  user,  the  explicit  boundary  condition  can  be  represented 
in  the  general  form: 
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2pi»rfl  “  Plidj 

Pt 

^Ph«u  ~  Pm 

Pn. 

^P^«wdl  P^M 

u 

u>dj 

V 

Vjjjj 

w 

T 

W«lj 

L 

be 

2T wall  ~  Tajj 

where 

Pw  *  Pm  ^  noncatalytic  wall 


s  calculated 


if  catalytic  wall 


(3-30) 


and 


Twtu  =  if  adiabatic  wall 

s  user-specified  constant  if  isothermal  wall 

The  catalytic  wall  boundary  conditions  are  discussed  below. 

Catalytic  Wall.  In  the  nonequilibrium  flow  over  a  body,  the  solid  surface  may  act  as  a  catalyst 
for  the  recombination  of  atoms  and  ions:  hence,  the  heat  transfer  at  the  surface  will  increase. 
Reentry  heating  data  derived  from  space  shuttle  missions  clearly  showed  the  significance  of 
nonequiUbrium  gas  chemistry  on  aerodynamic  heating  in  a  high-velocity,  low-density  flight 
legiine.  As  altimde  and/or  velocity  is  increased  above  shuttle  reentry  levels,  nonequilibrium 
effects  become  even  more  pronounced.  Consequently,  the  capability  to  model  surface  kinetics  is 
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essential  in  developing  a  flowfield  code  applicable  to  high-speed,  low-density  flows. 

The  HYUDA  code  models  wall  catalyticity  effects  as  part  of  the  solid  wall,  no-slip 
boundary  condition.  Wall  catalycity  is  represented  in  terms  of  species  catalytic  efficiencies, 
defined  as  the  ratio  of  the  number  of  species  recombining  at  the  wall  to  the  number  of  species 
arriving  at  the  wall.  Determination  of  the  species  catalytic  efficiency  is  dependent  upon 
experimental  measurements,  and  there  is  a  degree  of  uncertainty  as  to  the  best  possible  way  to 
represent  surface  catalytic  effects. 

The  HYUDA  code  implements  the  catalytic  boundary  condition  in  a  lagged  manner.  That 
is,  the  new  values  of  the  species  densities  at  the  wall  because  of  catalycity  are  calculated  at  the 
completion  of  each  computational  step,  rather  than  as  part  of  the  step.  The  governing  equation 
for  the  calculation  of  the  catalytic  wall  species  densities  is: 


P.Ys 

2-Ys 

V  J 

waU 

2n 

(3-31) 


where 

p  =  total  density, 

i  =  species  s  diffusion  velocity,  defined  in  eq.  (3-4), 

Ys  =  species  s  catalytic  efficiency, 

5,  =  -fl  for  polyatomic  species,  -  I  for  monatomic  species. 


with  R  being  the  Universal  Gas  Constant  and  M,  being  the  species  molecular  weight. 

Equation  (3-31)  was  derived  subject  to  the  following  conditions: 

•  No  blowing,  ablation,  or  porous  surface; 

•  Neglects  other  species  (i.e.,  assumes  zero  concentration  for  all  species  except  the  species 
being  calculated); 
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•  Uses  the  Lagmuir-Hinschelwood  recombination  mechanism; 

•  Assumes  first  order  surface  catalytic  recombination. 

Applying  the  definition  of  the  diffusion  velocity  at  the  wall,  the  left-hand-side  of  eq.  (3-31) 
becomes: 


r  -1 

■  „  M.  dX,' 

pqa 

L  J  wall 

P^»aa»  VI 

^mix  an 

W  J 

waU 


(3-32) 


where  is  the  mixture  diffusion  coefficient  for  species  s  (defined  in  section  4.0),  and 
M«i.  are  the  species  and  mixture  molecular  weights,  respectively,  and  dX,/dn  is  the  gradient  of 
the  species  mole  fraction,  with  X,  =  (Ps/M,)/(p/M,„ix)‘ 

Substituting  eq.  (3-32)  into  eq.  (3-31)  and  solving  for  dX^/dn  yields: 


'  ax,' 

-X, 

Ya5,  ' 

RfTwall 

.  an. 

wall 

^  J 

wall 

2-y, 

27C 

k  4 

(3-33) 


Expanding  the  gradient  term  on  the  left-hand-side  of  eq.  (3-33)  in  a  first  order  difference 
and  solving  for  the  molar  concentration  in  the  boundary  cell  yields: 


(3-34) 


where 


ax, 

dn 


J  waU 


is  given  by  eq.  (3-33). 
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Finally,  the  boundary  cell  species  density  is  obtained  from: 


P»hc  ~^*be 


Pbc 


^^mixbe 


M. 


(3-35) 


Turbulent  Wall.  Hypersonic,  high-altitude,  low-density,  low-Reynolds  number  flows  such  as 
those  of  interest  to  this  contract  are  primarily  laminar  flows.  In  the  late  phases  of  reentry, 
though,  transition  to  turbulence  may  occur.  The  turbulence  expected  for  the  cases  of  interest 
here  include  shear  layer  type  turbulence  on  the  windward  side  of  the  reentry  vehicle,  as  well  as 
vortex-separated  flows  on  the  leeside  of  the  vehicle. 

Windward  turbulence  of  the  shear  layer  type  is  expected  to  occur  as  a  shock/boundary  layer 
interaction,  such  as  might  be  found  at  a  deflected  control  surface  for  example.  This  shear  layer 
type  turbulence  is  relatively  simple  to  model,  even  with  basic  algebraic  turbulence  models  such 
as  the  Baldwin-Lx)max  model  (ref.  4). 

Leeside  turbulence  occurs  when  the  vehicle  angle  of  attack  is  sufflcient  to  separate  the  flow 
over  the  leeside  siufaces,  forming  vortices  that  can  strongly  influence  heating  patterns.  These 
leeside  flows  are  highly  sensitive  to  vehicle  geometry,  angle  of  attack,  and  Reynolds  number. 

The  Interim  Report  for  Phase  !  of  this  contract  (ref.  2)  discussed  the  turbulence  model 
selection  process  in  the  context  of  both  shear  flows  and  leeside  flows.  The  turbulence  model 
selected  was  the  Baldwin-Lomax  algebraic  turbulence  model  (ref.  4)  subject  to  the  modifications 
of  Degani-Schiff  (ref.  22),  and  also  including  a  wail  function  Qaw-of-the-wall)  formulation. 
Details  of  the  model  are  included  here. 

The  selection  of  the  Baldwin-Lomax  model  over  more  complex  models  was  in  part 
influenced  by  the  need  to  calculate  finite  rate  chemistry  with  its  inherently  large  number  of 
equations.  The  increased  accuracy  obtained  by  more  complex  models  was  not  considered 
sufficient  to  offset  the  increased  computing  time  required  for  higher-order  models,  especially 
when  considering  finite  rate  chemistry  calculations  with  large  numbers  of  species.  The  HYLDA 
code  is,  however,  configured  to  permit  incorporation  of  higher-order  turbulence  models. 

The  Baldwin-Lomax  model  calculates  turbulent  viscosity  according  to: 
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ify<ycK»t 


(3-36) 


where 


^turb  =pK^h2F^la)l  , 

^  J  inner 

I  )^nirbl  =  KCqjFwtkeFlCLEB  »  (3-37) 

^  J  outer 


Vaou  is  smallest  distance  normal  to  the  solid  wall  at  which  the  inner  and  outer 
turbulence  viscosities  are  equal. 

The  variables  in  eq.  (3-37)  are  defined  as  follows: 

p  «  density, 

K  =  von  Karman  constamt  (=0.41) 

h  =  distance  normal  the  wall 
(0  =  vorticity 
Cep  =  constant  (=1.6) 

and  the  Van  Ernest  (Fvd)  and  Klebanoff  (Fifi  ca)  functions  are  defined  as: 


PRLEB  = 


1+5.5 


CkleB  Ifi  I 


Ihl 


-1 


(3-38) 


The  variables  in  eq.  (3-38)  are  defined  as  follows: 


y+  = 


tttypwdi 

PwaU 
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Oj  =  friction  velocity  =  (Twau/pwau)'^  . 
Pwtii  =  wall  viscosity  , 

'Cwtu  =  wall  shear  stress  , 

A"^  =  constant  ( =  26 )  , 

Cifira  =  Klebanoff  constant  ( =  0.3  )  , 


h  =  the  value  of  h  at-j  I  h  I  I  co  I  Fyo 


flux 


The  wake  function  Fwake  is  the  minimum  of  two  functions: 


Fwake  ~  min^  Fwl»Fw2j 


(3-39) 


where 


I h I  leal Fvd 


max 


Fw2  = 


Cwkihlq: 


max 


1 

"  1  h  1 1  ct)  1  Fvd 

» 

J 

max 

(3-40) 


with  Cwk  =  0-25,  and  q^ax  equal  to  the  maximum  velocity  in  the  boundary  layer. 

The  Baldwin-Lomax  model  as  shown  in  eqs.  (3-36)-(3-40)  is  sufficient  for  the  wall- 
bounded  turbul-*nt  shear  layers  found  on  the  windward  side  of  hypersonic  vehicles.  As  discussed 
in  the  Phase  I  Interim  Technical  Report  (ref.  2),  however,  the  definition  of  the  length  scale  h  in 
eq.  (3-38)  leads  to  an  improper  length  scale  determination  when  calculating  leeside  flows.  To 
correct  the  deficiency,  the  Degani-Schiflf  modification  (ref.  22)  to  choose  the  first  maximum  of 

A 

the  function  determining  h  in  eq.  38)  has  been  incorporated. 

Wall  Functions  (Law-of-the-Wall).  Turbulence  length  scales  are  generally  orders  of 
magmtude  smaller  than  laminar  length  scales,  and  to  adequately  resolve  a  turbulent  boundary 
layer  (i.e.  grid  to  within  y*  <  5)  requires  a  large  number  of  grid  points,  and  a  large  amount  of 
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computer  time.  Wall  functions  such  as  the  law-of-the-wall  formulation  used  in  the  HYLDA 
code  are  used  to  relax  the  amount  of  grid  necessary  to  provide  the  correct  wall  value  of  viscosity. 
The  law-of-the-wall  model  predicts  the  wall  viscosity  using  empirical  correlations  from 
experinoental  studies  of  turbulent  boundary  layers.  These  correlations  are  applied  at  greater 
distances  from  the  wall  (typically  y'*’  ~  of  1(X)),  thereby  reducing  the  number  of  grid  points 
needed  to  adequately  resolve  the  turbulent  boundary  layer. 

The  governing  equation  is  based  on  the  following  expression  for  the  shear  velocity  (ux): 


ai  = 


In 


KU 


yPwall 


PwaU 


Ut 


(3-41) 


where 


Ux  =  shear  velocity  (defined  in  eq.(3-38))  , 

K  =  von  Karman  constant  (  =  0.41)  , 

C  =  5.1  . 

K  -  1.25k  =  0.5 125  for  flat  plates  , 

W  =  Cole's  tabulated  universal  wake  function  (taken  here  as  0)  , 
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J 


with 

_r  0.5(Y-l)Mi 
l-»-0.5(Y-l)Mi 

Equation  (3-41)  is  iterated  to  obtain  the  correct  shear  velocity  at  the  wall,  and  the  wall  viscosity 
is  then  extracted  from  this  value  as  follows; 
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Solid  Wall  (Free-slip),  Symmetry.  The  solid  wall  (free-slip)  and  symmetry  boundary 
conditions  ensure  that  all  fluxes  are  reflected  at  the  boundary; 


Vbc=RV^j 


(3-43) 


where  R  is  the  rotation  matrix,  defined  as: 
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with  /—  being  the  direction  cosines. 
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4.0  THERMOCHEMISTRY  MODELS 


The  HYLDA  code  contains  three  types  of  thermochemistry  models:  ideal  gas,  equilibrium 
air,  and  finite  rate  (nonequilibrium)  air.  In  general,  thermochemical  models  provide  the 
chemistry  source  terms,  thermodynamic  properties,  and  transport  properties  for  the  calculation 
scheme.  For  ideal  gases  and  equilibrium  air  the  source  terms  are  zero.  For  an  ideal  gas,  the 
thermodynamic  pressure  and  temperature  are  calculated  directly  firom  the  ideal  gas  relationships 
derived  from  eqs.  (3-S)  and  (3-8)  respectively.  The  equilibrium  air  model  is  empirical,  and 
provides  pressure,  temperature,  speed  of  sound,  and  ratio  of  specific  heats  ("f)  given  internal 
energy  and  density  using  Tannehill’s  real  gas  model  (ref.  5).  For  both  ideal  gases  and 
equilibrium  air,  the  transport  property  viscosity  is  obtained  using  Sutherland’s  Law,  and  the 
thermal  conductivity  is  determined  using  a  viscosity  and  Prandtl  number  relationship. 

For  finite  rate  chemistry  flows,  the  source  terms  are  nonzero,  the  thermodynamic  and 
transport  properties  depend  upon  the  properties  of  the  constituent  species.  Finite  rate  chemistry 
flows  require  detailed  models  to  account  for  nonequilibrium  effects  on  the  source  terms, 
thermodynamics,  and  transport  properties.  The  results  of  this  contract  effort  include  the 
development  of  an  air  chemistry  model,  as  well  as  a  transport  properties  model  and  species 
thermodynamics  model.  The  remainder  of  this  sections  details  these  models. 

4.1  Air  Chemistry  Model 

An  extensive  literature  search  was  performed  for  chemical  and  thermal  nonequilibrium 
experimental  data.  Several  sources  of  shock  tube  experiments  were  found  (refs.  23,  24)  which 
could  be  used  to  calibrate  air  nonequilibrium  chemistry  mechanisms.  The  data  were  considered 
sufficiently  accurate  to  compare  various  chemical  nonequilibrium  models  (i.e.,  one-temperature 
models).  However,  the  same  was  not  true  for  thermal  nonequilibrium  (multi-temperature) 
experimental  data  sets.  Very  few  data  sets  were  found  which  contained  useful  measurements  of 
thermal  nonequilibrium  processes.  Based  on  the  results  of  the  literature  search  we  decided  that 
only  one-temperature  air  chemistry  models  could  be  calibrated  to  engineering  accuracy  and 
would  be  considered  in  this  connect.  Thus,  rotation,  translation,  and  vibration  energy  modes 
were  assumed  to  be  in  equilibrium. 

Many  air  ionization  models  were  identified,  with  each  model  utilizing  different  chemical 
species,  reaction  paths,  rate  coefficients,  and  application  features.  For  instance,  a  model  may 
depend  on  one  (translation)  or  two  (translation  and  vibration)  temperatures,  may  have  up  to  IS 
chemical  species,  may  use  an  Arrhenius  or  non-Airhenius  reaction  rate,  and  may  have  as  many 
as  50-60  reaction  paths.  These  variations  make  a  thorough  analysis  extremely  tedious  and 
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difficult.  To  insure  a  feasible  analysis,  the  following  assumptions  and  restrictions  were 
employed; 

•  Only  Arrhenius  reaction  rates  were  considered. 

•  Chemical  reaction  rates  were  assumed  to  depend  only  on  the  translational  temperature  (T); 

i.e.,  one-temperature  reaction  rates. 

•  Chemical  reaction  rates  were  assumed  to  be  valid  up  to  17000  K. 

•  Argon  has  a  negligible  effect  in  high-temperature  air  reactions. 

A  variety  of  species  have  been  employed  in  these  air  chendstry  models;  however,  the 
species  that  can  be  utilized  in  any  calculation  is  also  dependent  on  the  available  species 
thermodynamic  and  transport  properties.  A  review  of  species  thermodynamic  and  transport 
properties  at  the  irtitiation  of  this  contract  concluded  that  at  least  rune  air  species 
(N2,C)2,N,0,N0^,N^,0'^,E"  )  could  be  modeled  These  species  appear  to  adequately  model 
high-temperature  air  in  most  engineering  problems;  hence,  the  thermochemistry  capabilities 
currency  included  in  HYLDA  employ  the  nine  species  mentioned  above.  HYLDA  can  easily 
acconunodate  additional  species  if  atomic/molecular  information  for  each  species  is  available. 

Of  the  many  models  examined,  five  air  ionization  models  utilized  an  Anfienius  reaction 
model  of  the  form  AiT^'expf-Ej/T)  where  Aj.Ni.Ei  are  rate  coefficients.  These  models  were 
developed  by  K.  Wray,  C.  Park  and  G.Menees,  S.  Kang  and  M.  Duim,  D.  Bittker,  and  M. 
Bortner  (refs.  6,  25-28).  The  Wray  model  is  well  known  and  represents  one  of  the  simplest  (18 
reaction  paths,  7  species)  ionization  models  conceived  The  Park/Menees  model  was  derived  to 
investigate  meteoroid  wakes  and  both  the  Kang/Durm  and  Bortner  noodels  were  employed  for 
flowfield  analyses  about  reentry  vehicles.  Finally,  the  Bittker  model  was  employed  to  analyze 
shock-tube  air  reactions.  Note  that  the  Wray  model  has  been  modified  from  its  original  form. 
Modifications  included  the  elimination  of  N2  +  O2  =  2  NO,  based  on  C^amac  and  Feinberg’s 
recommendation  (ref.  23)  and  the  assumption  that  Argon  has  a  negligible  effect  in  high 
temperature  air  reactions.  The  Park/Menees  model  (ref.  6)  has  been  abbreviated  to  include  only 
nine  air  species  and  is  the  finite-rate  chemistry  model  included  in  HYLDA.  The  Appendix  lists 
the  reaction  paths  and  rate  constants  for  the  Park/Menees  model.  For  details  of  the  chemistry 
OKxlel  selection  process,  refer  to  the  Pha.«e  I  Interim  Technical  Report  (ref.  2). 
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4J  Thermodynamic  Properties 

Species  thennodynamic  properties  such  as  entropy,  Gibbs  free  energy,  enthalpy,  and 
specific  heat  capacity  are  important  components  to  model  development  These  properties  must 
be  accurately  known  up  to  very  high  temperatures  (20000  K).  McBride  (ref.  8)  derived 
polynomial  relationships  for  the  thermodynamic  properties  of  the  nine  air  species  considered 
above.  McBride  utilized  the  following  relationships  for  the  specific  heat  capacity  (Cp),  enthalpy 
(H)  and  Gibbs  free  energy  (F): 

=  Zi  +  Z2T  +  ZaT^  +  Z4T^  +  ZsT*  (4-1) 

=  Zi  +  Z2T/2  +  ZaT^/S  +  Z4TV4  +  ZsT^/S  +  Z^IT  (4-2) 

K.  1 

-^  =  Zi  [1.0-loge(T)]  -  Z2T/2  -  ZaT^/fi  -  Z4TV12  -  ZsT^nO  +  ZsrT-Zn  (4-3) 

The  coefficients  for  eqs.  (4-1)  through  (4-3)  for  all  nine  air  species  can  be  found  in  the 

Appendix. 

4J  Transport  Properties 

The  transport  of  fluid  momentum,  thermal  energy,  and  mass  is  characterized  by  the 
coefficient  of  viscosity  (^),  theimal  conductivity  (X),  and  (concentration)  difrusion  D, 
respectively.  Viscosity  affects  the  growth  of  the  boundary  layer,  hence,  will  influence  the 
distance  (and  time)  particles  and  energy  must  travel  to  reach  the  wall  surface.  Conversely, 
thermal  conductivity  regulates  the  transfer  of  heat,  and  thus,  influences  particle  motion  (i.e., 
species  concentration  and  thermal  diffusion). 

All  nonequilibrium  chemistry  methods  require  the  computation  of  each  transport  property 
for  individual  species  in  the  reactive  flowfield.  A  foimula  is  then  used  to  obtain  the  mixture 
value  of  the  property.  This  procedure  is  depicted  in  figure  3,  which  also  shows  the  fimctional 
dependence  of  each  property  on  other  parameters. 
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Figure  3.  Transport  Properties  Methodology 

First-approximations  to  the  species  transport  properties  can  be  represented  as  follows: 


ilij 

(mi  +  mj)A[f 

(4-4) 

15  ^boltz 

■  4  A[?  ’ 

(4-5) 

_  hjjoitjT 

"  pa|/)  ’ 

(4-6) 

where  =  Cq 


2mimj 


JtkboitzT(mi  +  mj) 


V4 


(4-7) 


r  -  1  r  - 

c,  -  3  .  Q  -  3 


(4-8) 
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where  ^,A.,D  are  the  coefficients  of  viscosity,  thermal  conductivity  and  mass  diffusion;  m  is  the 
particle  mass,  kboitz  is  the  Boltzmann  constant,  P  is  the  pressure  and  is  the  collision 

cross-section  (aica)  nonnalized  with  respect  to  its  rigid-sphere  value.  Equations  (4-4)  through 
(4-8)  are  based  on  kinetic  theory  of  dilute  gases.  Complete  derivations  are  provided  in  ref.  29. 

All  species  transpon  properties  depend  on  the  collision  integral  cross-section,  defined  as  the 
probable  collision  area  between  species  s  and  t  The  collision  cross-section  is  dependent  on 
species  mass,  velocity,  temperature,  collision  deflection  angle,  and  a  potential  force  function. 

The  cross-sections  for  all  possible  collision  pairs,  among  the  nine  air  species,  have  been 
modeled  by  Nicolet  (ref.  7)  using  a  simple  logarithmic  polynomial  function  of  temperature.  The 
Appendix  contains  the  coefficients  for  eq.  (4-9)  for  all  non-coulomb  type  collisions  of  the  form: 

jtflij'"  =  Aij^loforr/1000]  -I-  B|j"'  (4-9) 

Coulomb  collisions  involving  two  charged  particles  have  been  modeled  by  Liboff  assuming 
an  unscreened  Coulomb  potential  with  Debye-length  cutoff: 

Once  the  individual  species  transport  properties  are  known,  the  mixture  value  for  each 
property  may  be  obtained  by  the  mixture  relationships  of  Brokaw  (ref.  29): 
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Finally,  the  reciprocal  of  the  mixture  diffusion  coefficient  is  approximately  equal  to  the 
mass  average  of  the  reciprocal  of  the  binary  diffusion  coefficients  (refs.  30,  31),  as  follows: 


i=N 

Z  XjMj 
i=l.i*j _ 

i=N  Xi 
Mmix  Z  '5” 


(4-15) 


The  individual  species  formulas  in  cqs.  (4-4)  through  (4-8),  collision  integral  cross-section 
expressions  in  cqs.  (4-9)  and  (4-10),  and  the  mixture  equations  (4-11)  through  (4-14)  represent 
the  relationships  used  by  HYLDA  to  determine  the  transport  properties. 
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5.0  CODE  DEMONSTRATION 


The  demonstration  cases  detailed  in  this  secdon  were  calculated  as  part  of  the  contract 
requirements  for  the  current  HYLDA  code  development  effort  These  cases  were  selected  by  the 
Air  Force  for  the  purpose  of  demonstrating  the  HYLDA  code’s  finite  rate  chemistry  capabilities, 
and  should  not  be  interpreted  as  an  attempt  to  fully  validate  the  HYLDA  code.  Full  validation  of 
the  HYLDA  code  will  require  many  calculations  over  a  wide  range  of  test  conditions,  a 
requirenoent  beyond  the  scope  of  the  current  contract  Ongoing  (ref.  32)  and  future  Air  Force 
work  will  address  the  HYLDA  code  validation  issue  in  detail. 

Demonstration  case  1  (section  S.l)  consists  of  a  calculation  of  the  flow  about  a  Mach  IS 
wedge-flat  plate  combination,  and  illustrates  the  HYLDA  code  capability  to  calculate  reacting 
air  flowfields.  Demonstration  case  2  (section  5.2)  calculates  the  reacting  air  flow  about  a  biconic 
geometry  at  Mach  6.9,  with  emphasis  on  ^^urface  heat  flux.  Demonstration  case  3  (section  S.3) 
calculates  the  reacting  flow  about  a  hemisphere,  with  emphasis  on  the  effects  of  wall  catalycity. 

While  these  demonstration  cases  presented  here  should  not  be  considered  full  validation 
cases,  they  do  serve  three  purposes:  (1)  to  denx}nstrate  the  HYLDA  code  finite  rate  chemistry 
capability,  (2)  to  assess  the  ability  of  the  HYLDA  code  to  model  wind  tunnel  tests,  and  (3)  to 
illustrate  the  need  for  new  hypersonic  wind  turmel  results  capable  of  being  used  to  validate  new 
codes  such  as  HYLDA. 


5.1  Casel:  Mach  IS  Wedge-Flat  Plate  Flow 

The  first  demonstration  case  was  based  on  the  wedge-flat  plate  experimental  studies  by 
Vidal  and  Stoddard  (ref.  33).  The  objective  of  the  experimental  studies  was  to  demonstrate  the 
magrtitude  of  nonequilibrium  effects  as  compared  to  viscous  effects.  The  geometry  selected  was 
a  flat  plate  at  zero  incidence  preceded  by  a  wedge  leading  edge.  Vidal  and  Stoddard  studied 
both  thin  (0.5  inches,  0.0127  meters)  and  thick  (4.0  inches,  0.1016  meters)  plates,  with  the  thin 
plate  intended  to  simulate  nonreacting  gas  flows,  and  the  thick  plate  providing  reacting  gas  data. 
Since  the  goal  of  the  HYLDA  code  is  to  calculate  reacting  gas  flowfields,  the  focus  of  this 
demonstration  case  was  on  the  thick  wedge  geometry. 

As  noted  in  the  introduction  to  this  section,  one  purpose  of  the  demonstration  cases  was  to 
illustrate  the  need  for  improved  wind  tunnel  tests  of  hypersonic  flowfields.  A  major  problem 
associated  with  wind  tunnel  tests  of  the  early  1960’s  was  the  difficulty  in  accurately  measuring 
and  quantifying  flowfield  and  surface  data.  These  experimental  difficulties  translate  into 
uncertainties  in  modeling  the  flows  numerically,  since  in  many  cases  the  input  parameters  for  the 
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computer  code  (including  HYLDA)  must  be  extracted  &om  given  experimental  data. 

Another  uncertainty  in  the  experimental  results  from  shock  tunnels  is  the  effect  that  the 
transient  response  of  the  tunnel  had  on  the  "steady  state"  .reacrng  flowfield.  The  ability  of  the 
HYLDA  code  to  accurately  model  shock  tunnel  data  depends  upon  the  influence  of  the  transient 
response  in  the  shock  tunnel  on  the  measured  "steady  state"  flowfield  conditions. 

The  results  presented  below  for  this  demonstration  case  assumed  the  wedge-flat  plate 
geometry  to  be  in  the  test  section  of  die  wind  tunnel  so  that  die  oncoming  flow  could  be 
represented  as  parallel  flow.  Comparison  of  the  numerical  and  experimental  results  showed 
discrepancies  in  the  pressure  field  about  the  geometry.  Invesdgadon  of  both  the  numerical 
results  and  the  experimental  results  revealed  that  the  test  geometry  was  not  located  in  the  test 
section,  but  rather  in  the  inlet  nozzle  portion  of  the  wind  turmel.  The  HYLDA  code  currently 
requires  a  uniform  freestream  initial  condition,  so  the  radial  flow  of  the  experimental  analysis 
could  not  be  replicated. 

In  light  of  the  discrepancy  between  model  locations  (and  thus  irutial  conditions), 
quantitative  comparison  of  the  numerical  and  experimental  results  is  not  strictly  possible. 
Qualitative  comparisons  of  the  results  are  possible,  though,  and  the  results  presented  here  show 
that  some  flow  features  are  relatively  insensitive  to  the  radial  component  of  velocity,  allowing  at 
least  a  qualitative  assessment  of  the  ability  of  the  HYLDA  code  to  model  the  reacting  gas  flow 
over  the  wedge-flat  plate  geometry. 

The  discussion  of  results  is  divided  into  three  parts:  geometry,  initial  conditions,  and 
results  and  discussion. 

Geometry.  The  geometry  used  for  this  case  consisted  of  a  4  inch  thick  (0.1016  meters)  flat 
plate  with  a  42.5-degree  half-angle  wedge  leading  edge.  The -geometry  is  shown  in  fig.  4: 


Figure  4.  Demonstration  Case  1  Geometry 
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Initial  Conditions.  Some  of  the  initial  conditions  required  by  HYLDA  (freestream  pressure 
and  temperature)  were  not  explicitly  provided  by  the  discussion  of  the  experimental  results  in 
ref.  33.  These  conditions  were  extracted  firom  the  data  provided,  as  follows: 


Given:  p/psL.  M«„  and  Re/inch; 

Extract:  P»,  T., 

•  Determine  the  fteestream  density  (p.,)  from  p/psL.  where  psL  is  the  sea  level  density. 


•  Determine  the  freestream  viscosity  (p..)  using  p...  U...  and  Re/inch. 


•  Iterate  for  the  freestream  temperature  (T*)  using  Sutherland’s  Law  and  the  p,.  calculated 
above: 


P-(T«)  = 


1.458(10~^)T1:^ 
T«+ 110.3 


•  Calculate  the  freestream  pressure  (P«)  using  the  freestream  values  of  density  and 
temperature  and  the  mixture  gas  constant  (for  the  air  chemistry  case  here,  the  mixture  gas 
constant  was  287.88  J/(kg-K)). 


The  results  below  are  presented  for  two  of  the  large-scale  wedge  cases  studied  by  Vidal  and 
Stoddard,  namely,  the  low  stagnation  temperature  case  and  the  high  stagnation  temperature  case 
1  tests.  With  the  procedure  outlined  above,  and  using  the  appropriate  parts  of  table  1  from  ref. 
33,  the  initial  conditions  used  here  are  as  follows: 

Tabk  5.  DcBonstration  Case  1  Initial  ConditioBS 


HMovair  oondiliaat* 

Amhienl  n  model  location" 

Cua 

T.dUC) 

P,  (pdaJfAn^) 

P-A>a. 

U.  (ft/ijnA) 

Macfa# 

Ra/InJteAm 

P.  (pdaJIAi^) 

T.nuo 

LowMpL 

3922.3390 

3820,  Z6(I0’) 

6.01(l(r*) 

14900,4343 

13.0 

9210,363600 

2.1(10^x  M-43 

322,179 

Higbianp. 

16700, 8J(10’) 

2J8(l(r*) 

9370,2836 

16.3 

2680, 103300 

16(l<r*),  18 

133,83 

rikahwiil  frooi  Miic  ooadiiiaM  u  modd  locaiioa,  d«a  fiam  uUe  1  (nf.  33) 
**  ealeaiHad<foairiMmircaadiliaaf,dMafranitaUi  I  (icf.33) 

***  ciknliiiil  from  nf.  33  Miibieat  oondiiiaai 


Additionally,  the  geometry  wall  temperature  was  assumed  to  be  equal  to  the  ambient  pre-test 
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temperature  of  300  K  (540  R). 


R<^Its  and  Discussion.  Using  the  initial  conditions  in  table  5,  the  HYLDA  code  was  used  to 
calculate  the  flowfield  about  the  wedge-flat  plate  geometry  for  both  the  low  stagnation 
temperature  case  and  the  high  stagnation  temperature  case  1  from  the  Vidal  and  Stoddard  report 
(ref.  33).  The  results  and  discussions  for  these  cases  will  be  described  separately. 

Low  Stagnation  Temperature  Results 

The  grid  and  the  computational  domain  for  the  low  stagnation  temperature  case  are  shown 
in  fig.  5.  The  grid  consists  of  80  cells  in  the  streamwise  direction,  50  cells  in  the  radial  direction, 
and  1  cell  in  the  spanwise  direction.  Only  one  point  is  required  in  the  third  direction  since  the 
flow  is  essentially  two-dimensional.  The  grid  is  clustered  near  the  nose  of  the  wedge  to  resolve 
the  shock,  and  also  clustered  near  the  body  to  resolve  the  viscous  layers. 

The  low  stagnation  temperature  calcularion  was  accomplished  in  112  steps  with  a 
maximum  CFL  number  of  500.  Figs.  6  through  8  show  the  Mach,  temperature,  and  y  contours 
for  the  flowfield,  where  y  is  the  ratio  of  specific  heats.  The  Mach  and  temperattirc  contours 
qualitatively  show  the  expected  flow  features,  including  the  shock  attached  to  the  leading  edge 
of  the  wedge  and  the  temperature  rise  through  the  shock  and  subsequent  cooling  to  the  wall 
temperature.  The  y  contours  are  useful  to  show  where  finite  rate  chemistry  effects  are  largest. 
As  expected,  the  greatest  degree  of  reacting  chemistry  occurs  in  the  high  temperature  region 
behind  the  shock  along  the  wedge  face. 

Vidal  and  Stoddard  (ref.  33)  included  a  Schlieren  photograph  for  this  low  temperature  case. 
Although  the  correlation  is  not  exact,  density  profiles  from  the  numerical  calculation  can  be  used 
to  compare  to  the  Schlieren  photograph.  The  comparison  of  numerically-calculated  density 
contours  and  experimentally-produced  Schlieren  photograph  is  shown  in  fig.  9(a,b),  with  part  (a) 
being  the  Schlieren  photograph,  and  part  (b)  being  the  numerical  results.  The  figures  show  that 
the  HYLDA  code  predicts  the  shock  wave  angles  on  both  the  upper  and  lower  surface.  The 
boundary  layer  on  the  flat  plate  is  more  easily  observed  in  figs.  6  through  8. 

As  noted  in  the  introduction  to  this  section,  direct  comparison  of  quantitative  results  is  not 
appropriate  because  of  the  difference  in  initial  conditions  for  the  experiment  and  the  calculation. 
The  comparison  of  surface  pressures  shown  in  fig.  10,  though,  demonstrates  that  the  pressure 
field  is  not  greatly  effected  by  the  radial  component  of  velocity  in  the  experimental  study. 
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Figure  5.  Demonstration  Case  1  Mesh  for  Low  Stagnation  Temperature  Case 
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Figure  6.  Demonstration  Case  1  Lour  Stagnation  Temperature  Case  •  Macta  Contours 
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Figure  7.  Demonstration  Case  1  Low  Stagnation  Temperature  Case  -  Temperature  Contours 
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Figure  8.  Demonstration  Case  1  Low  Stagnation  Temperature  Case  -  y  Contours 
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Figure  9a.  Demonstration  Case  1  Low  Stagnation  Temperature  Case  •  Schlieren  (ref.  33) 


51 


Figure  9b.  Demonstration  Case  1  Low  Stagnation  Temperature  Case  •  Density  Profiles  for 
Comparison  With  Schlieren  (fig.  9a). 
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Figure  10.  Demonstration  Case  1  Low  Stagnation  Temperature  Case  -  Pressure 
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High  Stagnation  Temperature  Results 

The  grid  and  the  computational  domain  for  the  high  stagnation  temperature  case,  shown  in 
fig.  11,  are  similar  to  those  used  for  the  low  stagnation  temperature  case.  Using  the  shock 
structure  from  the  low  stagnation  temperature  case  as  an  example,  the  grid  for  the  high 
sts^gnation  temperature  case  was  reduced  from  the  80x50x1  used  for  the  low  temperature  case  to 
a  70x35x1  grid,  with  grid  points  removed  primarily  from  the  freestieam  section  of  the  original 
80x50x1  grid.  Again,  only  one  cell  is  required  in  the  third  direction  since  the  flow  is  essentially 
two-dimensional.  The  grid  is  still  clustered  near  the  nose  of  the  wedge  to  resolve  the  shock,  and 
also  clustered  near  the  body  to  resolve  the  viscous  layers. 

The  high  stagnation  temperature  calculation  was  accomplished  in  448  steps  with  a 
maximum  CFL  number  of  1000.  Figs.  12  through  14  show  the  Mach,  temperature,  and  y 
contours  for  this  case,  where  y  is  the  ratio  of  specific  heats.  In  comparing  to  the  low  stagnation 
temperature  case  discussed  above  (figs.  6  through  8),  we  see  that  the  features  of  the  flowfield  are 
qualitatively  similar.  The  shock  is  closer  to  the  body  in  the  high  stagnation  temperature  because 
the  Mach  number  is  slightly  higher.  The  maximum  temperature  for  this  case  is  near  5200K 
compared  to  23(X)K  for  the  low  stagnation  temperature  case. 

The  y  contours  provide  an  indication  of  the  nonequilibrium  effects  of  this  higher 
temperature,  higher  Mach  number  test  case.  In  comparing  to  the  corresponding  low-temperamre 
y  contours,  we  see  that  the  nonequilibrium  flow  in  the  low  stagnation  temperature  case  leads  to  a 
y  of  1.289,  while  the  higher  temperature  casc  leads  to  a  y  of  1.270,  showing  that  the 
nonequilibrium  effects  are  slightly  larger  in  the  high  stagnation  temperature  case.  In  both  cases, 
the  maximum  nonequilibrium  effects  occur  behind  the  shock  along  the  wedge  face. 

As  noted  earlier,  direct  comparisons  are  not  strictly  applicable  here  because  of  the 
difference  in  iiutial  conditions  between  the  experiment  and  the  calculation.  However,  the 
difference  in  conditions  is  not  large  enough  to  completely  preclude  at  least  a  qualitative 
comparison  of  the  surface  data  obtained  from  these  tests.  Figs.  15  and  16  show  a  comparison  of 
the  experimental  surface  pressure  and  heat  flux  measurements  with  the  HYLDA  results  for  those 
same  quantities.  Also  included  in  the  figures  are  the  results  of  an  ideal  gas  calculation  of  the 
high-temperature  case  using  HYLDA. 

Fig.  15  shows  the  surface  pressures  from  the  experimental  study  as  well  as  the  numerical 
results  from  the  ideal  gas  and  finite  rate  calculations.  The  figure  shows  that  the  largest 
discrepancy  between  the  experimental  results  and  the  calculated  results  occurs  at  the  single  data 
point  on  the  wedge  face  (x/t^edge  ”  0-^)»  with  the  ideal  gas  value  falling  between  the 
experimental  data  point  and  the  finite  rate  chemistry  calculation.  Less  discrepancy  is  found 
along  the  flat  plate  afterbody,  but  again  the  ideal  gas  calculation  falls  between  the  experimental 
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and  finite  rate  chemistiy  calculation  results. 

Fig.  16  shows  the  surface  heat  flux  along  the  body,  again  for  the  experimental  results  and 
both  the  ideal  gas  and  finite  rate  chemistry  calculations.  This  figure  shows  greater  sensitivity  to 
the  difference  in  fieestream  conditions,  with  three  distinct  and  widely  varying  heat  flux  results. 
The  ideal  gas  results  show  the  greatest  surface  heat  flux,  as  expected.  The  finite  rate  chemistry 
calculation  shows  the  correct  trend  of  decreasing  heat  flux  since  the  nonequilibiium  chemistry 
effects  extract  some  of  the  energy  from  the  flow.  The  experimental  investigation,  however, 
shows  a  significantly  lower  degree  of  heat  flux. 

Possible  explanations  for  the  wide  differences  in  heat  fluxes  shown  in  fig.  16  include  the 
differences  in  the  freestream  velocity  field,  the  transient  response  of  the  heat  flux  gauges  used  in 
the  experiment,  and  lack  of  resolution  of  the  thermal  boundary  layer  in  the  numerical 
calculations.  Additional  study  is  required  to  determine  the  exact  nature  of  the  discrepancy. 
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Figure  11.  Demonstration  Case  1  Mesfa  for  High  Stagnation  Temperature  Case 
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Figure  12.  Demonstratiou  Case  1  High  Stagnation  Temperature  Case  •  Mach  Contours 
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Figure  13.  Demonstration  Case  1  High  Stagnation  Temperature  Case  •  Temperature  Contours 
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Figure  15.  Demonstration  Case  1  High  Stagnation  Temperature  Case  -  Pressure 
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Figure  16.  Demonstration  Case  1  High  Stagnation  Temperature  Case  •  Heat  Flux 


S2  Case  2:  Mach  6.9  BIconic  Flow 


Validation  case  2  investigated  the  surface  heat  flux  along  a  straight  biconic  at  Mach  6.9. 
The  biconic  geometry  is  depicted  in  fig.  17  and  was  oriented  at  an  angle  of  attack  of  zero  degrees 
to  the  flow.  The  fluid  medium  was  air  and  the  experiment  was  performed  at  NASA  Langley’s 
Expansion  Tube  Facility  (ref.  34).  The  frtestream  velocity,  pressure,  and  temperature  are 
shown  in  Table  6  below.  Nonequilibrium  effects  near  the  stagnation  region,  if  any,  are  expected 
to  be  small  and  a  seven-species  finite-rate  air  chemistry  model  (ref.  35)  was  employed  to 
evaluate  HYLDA’s  nonequilibrium  heat  transfer  capabilities. 

Table  6.  Demonstration  Case  2  Initial  Conditions 


u. 

P. 

T. 

Fluid 

5326  m/i 

2182  NW 

1604K 

Air 

The  computational  grid  for  this  test  case  is  shown  in  fig.  18.  The  grid  included  219  cells 
along  the  wetted  length,  30  cells  normal  to  the  body  and  1  cell  in  the  circumferential  direction. 
A  high  grid  density  was  employed  along  the  body  and  near  the  wall  since  the  primary  interest  of 
this  case  was  to  accurately  predict  the  surface  heat  flux.  Numerical  simulation  of  the  expansion 
tube  test  was  performed  with  the  measured  wall  temperature  distribution  shown  in  fig.  19.  This 
distribution  was  input  into  the  HYLDA  code  by  modifying  the  face  5  specified  temperature  wall 
boundary  condition  through  the  use  of  a  data  statement  This  modification  was  specific  to  this 
particular  demonstration  case  and  was  not  delivered  with  the  HYLDA  code. 

Fig.  20  shows  a  comparison  of  the  predicted  wall  heat  flux  with  measured  values  at  the  0- 
and  180-degree  meridians.  The  predictions  are  generally  lower  than  measurements  along  both 
conic  sections.  A  parabolized  Navier-Stokes  prediction  (ref.  34)  yielded  slightly  improved 
agreement  with  measurements  along  the  first  cone;  however,  very  similar  results  were  obtained 
along  the  second  conic  section.  Ref.  34  suggests  several  experimental  uncertainties  that  also 
may  have  influenced  the  experimental  accuracy  such  as  uncertainties  in  the  freestream 
conditions  and  the  establishment  of  a  non-uniform  flowfield.  Computational  uncertainties 
regarding  grid  density  effects  must  also  be  considered.  Additional  study  is  required  to  determine 
the  exact  nature  of  the  anomalies  in  the  predictions. 
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q  (MW/m^) 


Axial  Distance  (x/L) 


Figure  20.  Heat  Flux  Comparison  for  the  Straight  Biconic 
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5  J  Case  3:  Mach  5  Catalytic  Blunt  Sphere  Flow 

The  third  demonstration  case  was  based  on  an  arc  jet  study  of  a  catalytic  hemisphere  in  a 
Mach  5  reacting  Nitrogen  flow  (ref.  36).  The  purpose  of  this  case  was  to  compare  the  stagnation 
point  heat  flux  calculated  by  HYLDA  with  the  experimentally-obtained  value.  This  single  point 
of  comparison  proved  difficult  to  match  numerically,  so  the  study  of  this  case  was  expanded 
beyond  the  single  data  point  comparison  to  include: 

•  Comparison  of  stagnation  point  heat  flux, 

•  Effects  of  wall  catalycity  on  stagnation  point  heat  flux, 

•  Effects  of  grid  clustering  on  stagnation  point  heat  flux,  and 

•  Mixture  effects  on  stagnation  point  heat  flux. 

This  expanded  approach  also  provides  insight  into  the  code  sensitivities  to  both  grid  and  physical 
effects. 

Geometry  and  Initial  Conditions.  The  geometry  used  for  this  case  was  a  5.08-cm  (2.0- 
inch)  nose  radius  hemisphere.  The  initial  conditions  provided  by  tables  1  and  2  in  ref.  36, 
supplemented  by  additional  parameters  provided  upon  request  by  the  ref.  36  authors  were  as 
follows: 


Table  7.  Demonstration  Case  3  Initial  Conditions 


M.=5.4 

U.  =  5797m/«ec  (19020  ft/jec) 
P.  =  l.aUKT*)*!!!! 

T_  =2920  K  (5256  R) 

Xn.  =  0.6 
X„=0,4 


As  was  the  case  for  demonstration  cases  1  and  2,  it  was  difficult  to  determine  the  required 
HYLDA  input  parameters  firom  the  discussion  of  the  experiment  in  ref.  36;  thus,  there  is  some 
uncertainty  as  to  the  accuracy  with  which  the  numerical  calculation  models  the  actual 
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experiment.  This  difficulty  to  ensure  that  the  numerical  calculation  accurately  represents  the 
experimental  conditions  is  a  problem  that  should  be  addressed  when  performing  any  new 
hypersonic  wind  tunnel  tests  so  as  to  provide  adequate  infonnation  about  the  test  to  allow 
accurate  numerical  modeling. 

The  grid  consisted  of  a  54x30x1  mesh,  with  54  cells  in  the  streamwise  direction,  30  cells  in 
the  radial  direcdon,  and  1  cell  in  the  circumferential  direction  since  the  flow  was  essentially 
axisymnaetric.  The  grid  is  shown  in  fig.  21  below: 


Figure  21.  Demonstration  Case  3  Mesh 

Results.  As  noted  above,  the  actual  comparison  to  experimental  data  was  limited  to  a 
single  data  point  The  poor  comparison  with  the  experimental  results  resulted  in  an  expansion  of 
the  case  to  include  studies  on  the  effects  of  grid  and  physical  parameters  on  the  heat  flux  results. 
The  four  parts  to  the  expanded  study  are  discussed  separately  below. 

Comparison  of  Stagnation  Point  Heat  Flux.  The  first  attempt  to  compare  with  the 
experimental  data  consisted  of  a  single  data  point  for  the  stagnation  point  heat  flux  on  the 
hemisphere,  using  a  wall  catalytic  efficiency  of  0.01.  The  experimental  data  point  was  obtained 
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from  fig.  2  of  ref.  36,  and  the  numerical  data  point  was  obtained  using  HYLDA.  For  the 
specified  catalytic  efficiency  of  0.01,  the  experimental  value  of  the  stagnation  point  heat  flux  was 
0.49(10*)  W/m^,  and  the  calculated  value  was  1.77(10*). 

There  is  a  large  difference  between  the  experimental  and  numerical  values,  and  it  is  not 
known  if  the  discrepancy  is  due  to  an  inaccurate  modeling  of  the  flow  conditions,  or  an  error  in 
the  calculation.  To  study  the  possibilities  for  error,  addidonal  data  points  were  considered.  In 
this  case,  ref.  36  provided  some  indication  of  the  stagnation  point  heat  transfer  expected  for 
various  other  catalytic  efficiencies.  These  data  points  were  considered  next 

The  same  geometry  and  grid  were  used  to  calculate  the  stagnation  point  heat  flux  for 
catalytic  efficiencies  of  y*  =  0.00,  0.01,  0.50,  and  1.00.  The  results  of  these  calculations  are 
shown  in  fig.  22  below,  with  corresponding  experimental  data  points  obtained  from  fig.  2  in  ref. 
36. 


Catalytic  Efficiency,  y. 


Figure  22.  Stagnation  Point  Heat  Flux  for  a  Range  of  Catalytic  Efficiencies 

Fig.  22  shows  that  the  HYLDA-predicted  stagnation  heat  flux  compares  more  closely  to  the 
experimental  values  at  the  higher-catalyticity  values,  while  overpredicting  considerably  the 
low-catalyticity  cases.  The  dotted  lines  in  the  figure  were  added  only  to  distinguish  the  data 
points.  Because  the  noncatalytic  case  is  the  most  basic  case,  and  because  the  largest 


discTq)ancies  occuired  for  the  noncatalytic  case,  the  noncatalytic  case  was  chosen  for  further 
study. 

Effects  of  Grid  Oustering  on  Heat  Flux.  The  third  test  consisted  of  calculating  the 
stagnation  point  heat  flux  for  the  noncatalytic  case  using  grids  with  varying  clustering  of  grid 
points  near  the  wall  in  the  stagnation  region.  There  is  strong  evidence  to  suggest  that  grid 
clustering  near  the  wall  can  significantly  effect  the  calculated  heat  flux  (ref.  37).  Fig.  23  (a-c) 
shows  the  nose  region  grids  for  the  three  different  grid  clusterings  used  for  this  test,  ranging 
firom  coarse  (grid  a)  to  fine  (grid  c). 


a.  Coarse  b.  Medium  c.  Fme 

Figure  23.  Grid  Clustering  Schemes 

Table  8  shows  the  results  of  the  analysis,  including  the  grid  spacing  of  the  first  cell  at  the 
stagnation  point  and  the  associated  calculated  heat  flux  at  the  stagnation  point  Also  included  in 
the  table  is  the  noncatalytic  value  firom  fig.  2  in  ref.  36. 

Table  8.  Grid  Clustering  Effects  on  Stagnation  Point  Heat  Transfer 


Grid 

Fim  CeU  Height  (m) 

Sugnadon  Hen  Bux  (W/m^) 

• 

6i21(Kr‘) 

1.37(10*) 

b 

171(1 0-‘) 

1.44(10*) 

c 

1.43(!(r’) 

1.61(10*) 

Expeiimeiit 

- 

0J0(10*) 
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Table  8  shows  that  grid  spacing  near  the  wall  does  have  some  effect  on  the  stagnation  point 
heat  flux,  but  in  this  case  it  is  not  significant  enough  to  be  considered  the  only  cause  of  the 
discrepancy  between  the  calculated  and  experimental  results. 

In  a  final  attempt  to  quantify  the  difference  between  experimental  and  numerical  results,  the 
effect  of  the  firestream  mixture  properties  was  investigated. 

Mixture  Effects  on  Heat  Flux.  Since  it  is  often  difficult  to  deteraune  accurately  the 
freestream  conditions  in  hypersonic  wind  tunnel  tests,  one  possible  source  of  inaccuracy  could 
be  related  to  the  presence  of  a  different  fireestream  than  that  reported  for  the  wind  mnnel  test. 
This  final  calculation  was  directed  at  studying  the  effect  that  a  different  freestream  mixture 
would  have  on  the  stagnation  point  heat  flux.  This  calculation  used  the  same  54x30x1  grid 
shown  in  fig.  21,  and  the  same  conditions  shown  in  table  7,  with  the  exception  that  the 
freestream  was  assumed  to  consist  of  pure  N2  rather  than  the  mixture  of  N2  and  N  as  reported  in 
ref.  36.  The  results  for  this  calculation  are  plotted  in  fig.  24; 
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Figure  24.  Freestream  Effects  on  Stagnation  Point  Heat  Flux 
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Fig.  24  is  the  same  as  fig.  22,  with  the  new  data  points  for  the  pure  N2  freestream  added. 
Again,  dotted  lines  were  added  to  the  data  points  from  fig.  22  to  distinguish  data  points.  Dotted 
lines  were  not  drawn  for  the  new  data  points,  since  only  two  catalytic  efficiencies  (0.00  - 
noncatalytic,  and  l.(X)  -  fuUy-catalytic)  were  calculated.  Fig.  24  shows  that  the  freestream 
effects  for  this  calculation  were  not  significant  at  the  fi^y-catalytic  value,  decreasing  slightly  the 
previous  HYLDA-calculated  value  using  the  ref.  36  reported  fireestream.  The  new  value  for  the 
noncatalytic  case  increased  over  the  previous  HYLDA  calculation,  leading  to  worse  comparison 
with  the  experimental  value.  In  general,  though,  the  effects  of  the  different  freestream  were  not 
large,  and  cannot  be  interpreted  as  the  cause  of  the  discrepancy  between  the  numerical  and 
experimental  results. 

The  goal  of  this  demonstration  case  was  initially  one  of  comparing  stagnation  point  heat 
flux  for  a  catalytic  blunt  sphere.  However,  the  discrepancy  between  the  numerical  and 
experimental  results  led  to  a  study  of  various  effects  on  heat  flux,  including  catalyticity  effects, 
grid  effects,  and  fieestreara  mixture  properties  effects.  While  none  of  these  effects  is  the  cause 
of  the  discrepancy,  each  of  these  studies  did  effect  the  solution  to  some  degree.  Further,  the 
uncertainty  in  the  exact  degree  to  which  the  numerical  calculation  accurately  represented  the 
experimental  test  is  not  clear,  and  precludes  additional  testing  on  this  particular  case. 
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6.0  CONCLUSIONS  AND  RECOMMENDATIONS 


The  HYLDA  code  developed  for  this  contract  is  capable  of  predicting  trends  in  hypersonic 
flowfields  including  real  gas  effects  on  heating  rates,  pressures,  forces,  and  moments  on  complex 
shapes  in  low-density  flows.  The  code  is  based  on  a  fully-coupled  Navier-Stokes/finite  rate 
chemistry  algorithm,  using  an  implicit  solution  algorithm  and  a  one  temperature  chemistry 
model.  The  species  thermodynamics  properties  model,  and  the  transport  properties  model  have 
been  tailored  to  fit  the  requirements  for  calculating  finite  rate  air  chemistry  up  to  Mach  25. 

The  HYLDA  code  is  a  versatile  code  capable  of  calculating  ideal  gas  flows,  equilibrium  air 
flows,  and  finite  rate  reacting  air  flows.  The  code  is  written  in  a  top-down  format  with  a  flexible 
data  structure,  both  of  which  facilitate  the  incorporation  of  new  technology  into  the  code  as  it 
becomes  available.  The  chemistry  model,  species  thermodynamics  model,  and  the  transport 
properties  model  are  also  flexible,  allowing  different  chemistry  models  to  be  installed  with 
minimal  changes. 

The  demonstration  cases  included  in  this  report  are  not  an  attempt  to  fully  validate  the 
HYI  DA  code.  They  do,  however,  demonstrate  some  of  the  capabilities  of  the  code  to  calculate 
finite  rate  chemistry  flows,  while  at  the  same  time  accentuating  the  need  for  new  wind  tunnel 
data  on  hypersonic  flows  to  support  the  validation  efforts  of  new  Navier-Stokes/finite  rate 
chemistry  codes.  The  demonstration  cases  provide  insight  into  the  approach  used  to  calculate 
hypersonic  reacting  flows,  and  provide  details  of  the  difficulty  in  trying  to  numerically  reproduce 
experimental  results.  Full  validation  of  the  HYLDA  code  will  require  calculation  of  a  large 
number  of  cases  over  a  wide  range  of  conditions,  and  will  require  a  significant  amount  of  time  to 
complete.  The  Air  Force  has  already  begun  efforts  to  validate  the  HYLDA  code. 

Recommendations  for  future  work  with  the  HYLDA  code  can  be  grouped  into  four 
categories: 

•  Validation  studies, 

•  Boundary  condition  extensions, 

•  Algorithm  extensions,  and 

•  Alternate  and  improved  chemistry  models. 

Validation  studies  should  include  both  numerical-to-experimental  comparisons  and 
numerical-to-numerical  comparisons.  Meaningful  numerical-to-experimental  comparisons  wiU 
require  careful  analysis  of  existing  experimental  hypersonic  analyses,  and  greater  emphasis  on 
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new  and  improved  methods  for  obtaining  accurate  hypersonic  experimental  data  using 
techniques  such  as  those  described  in  ref.  1.  Numerical- to-numerical  comparisons  will  provide  a 
measure  of  the  effects  of  different  numerical  techniques,  in  particular  the  chemistry  models 
being  used.  As  new  models  become  available,  codes  such  as  HYLDA  that  are  able  to  rapidly 
and  efficiently  incorporate  new  models  will  facilitate  the  analysis  of  the  models. 

Boundary  condition  extensions  for  the  HYLDA  code  should  include  the  addition  of 
multizone  capability  to  allow  calculation  of  more  complex  geometries  than  those  being 
considered  in  the  current  contract  effort.  Extensions  to  the  physical  boundary  condition  models 
should  include  porous  walls  (to  include  transpiration  and/or  suction)  and  improved  catalytic 
models.  The  flexible  data  structure  of  the  HYLDA  code  will  enable  incorporation  of  these  and 
other  boundary  condition  improvements  without  requiring  major  code  revisions. 

Algorithm  improvements  to  HYLDA  should  include  incorporation  of  new  developments 
intended  to  make  the  algorithm  more  stable  (to  save  computer  time)  and  to  increase  code 
robusmess  (to  extend  the  range  of  applicability).  These  improvements  should  consider 
alternative  methods  of  flux  splitting,  improved-  use  of  multiprocessing  capabilities,  and 
adaptive-chemistry  capability.  Again,  the  HYLDA  data  structure  should  facilitate  the 
incorporation  of  new  technology. 

Chemistry  models  are  the  most  dynamic  models  in  use  in  the  HYLDA  code.  Many 
special-purpose  chemistry  models  exist,  and  could  be  incorporated  into  HYLDA.  Many  other  air 
chemistry  models  exist,  and  these  models  could  be  incorporated  into  HYLDA  to  facilitate  direct 
comparisons  between  HYLDA  and  other  finite  rate  air  chemistry  codes.  Of  particular  interest 
should  be  additional  transport  properties  models  and  species  thermodynamics  databases. 

The  flexibility  of  the  HYLDA  Navier-Stokes/finite  rate  chemistry  code  should  make  it  a 
valuable  design  tool  for  the  next- generation  of  hypersonic  vehicles  being  developed. 
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ABBREVIATIONS 


AFWL 

Air  Force  Wright  Laboratories 

CPU 

central  processing  unit 

fig 

figure 

G-S 

Gauss-Seidel 

HYLDA 

Hypersonic  Low  Density  Analysis  code 

ref 

reference 

3D 

three-dimensional 
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NOMENCLATURE 


English 


_n+l  _n+l 

ai  ,82  ,83 


84,85,8? 


A,  B,  C,  D 


A^.B^.C^ 


A-.  B",  cr 


Ckleb 


implicit  coefficient  matrices,  eq.  (3-18) 
explicit  coefficient  matrices,  eq.  (3-18) 
inviscid  flux  Jacobians,  eq.  (3-13) 

positive  eigenvalue  flux  Jacobians,  eqs.  (3-14),  (3-16),  (3-17) 

negative  eigenvalue  flux  Jacobians,  eqs.  (3-14),  (3-16),  (3-17) 

constant,  eq.  (3-38) 

constant,  eq.  (3-41) 

constants  defined  by  eq.  (4-8),  i=l,2 

constant,  eq.  (3-37) 

Klebanoff  constant,  eq.  (3-38) 
specific  heat  at  constant  pressure 
specific  heat  at  constant  volume 
constant,  eq.  (3-40) 

species  s  mixture  diffusion  coefficient,  eq.  (3-4) 

“7  'I  ”7 

total  energy  per  unit  volume  =  p  +  V2  u  -(-  v"^  -1-  w  ,  eq. 
(3-2) 
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Ci 


Ci 


S 


®>elec 


Ci 


rot 


■  trans 


®‘vib 

f,(T) 

F 


r.G'.H' 
F,  G,H 


Fi.Gi.Hi 
Fy,  Gy,  Hy 
F^,  Gt^, 


Fkleb 


Fvd 


F 


wake 


wl 


total  internal  energy  per  unit  mass,  eq.  (3-6) 

species  s  internal  energy,  eq.  (3-7) 

electronic  internal  energy,  eq.  (3-6) 

heat  of  rormation,  eq.  (3-6) 

rotational  internal  energy,  eq.  (3-6) 

translational  internal  energy,  eq.  (3-6) 

vibrational  internal  energy,  eq.  (3-6) 

curve-fit  function  of  specific  heat  at  constant  pressure 

Gibbs  free  energy,  eq.  (4-3) 

X-,  y-,  z-diiection  flux  vectors,  eqs.  (3-1),  (3-2) 

^-direction  (transformed)  flux  vectors,  eqs.  (3-9),  (3-10) 
inviscid  flux  vectors,  eq.  (3-11) 
viscous  flux  vectors,  eq.  (3-11) 
transformed  viscous  flux  vectors,  eq.  (3-17) 

Klebanoff  intermittency  function,  eq.  (3-38) 
van  Driest  function,  eq.  (3-38) 
wake  function,  eq.  (3-39) 
wake  function  1,  eq.  (3-38) 
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wake  function  2,  eq.  (3-38) 


distance  nonnal  to  the  wall,  eq.  (3-37) 
mixture  heat  of  fonnation 
species  s  heat  of  formation 
species  s  enthalpy 
enthalpy,  eq.  (4-2) 
identity  matrix,  eq.  (3-13) 

Jacobian 

Boltzmann  constant,  eq.  (4-6) 
direction  cosines  (m=l-3,  j=l-3) 
mass  of  particle  i,  eq.  (4-6) 
species  s  molecular  weight 
Mach  number 

time  level 

electron  number  density,  eq.  (4-10) 
normal  direction 

matrix  relating  nonconservative  and  conservative  solution  vecton, 
cq. (3-13) 
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N, 

P>P 

q 

q 


qe 


qmu 


R 


A 

R 

Rs 


Re 

t 


^wedge 


T 


(u,v,w) 


u 


number  of  species 
pressure 

<u,v,w>,  eq.  (3-3) 
heat  flux 

edge  velocity,  eq.  (3-41) 

<  u,  ,  V, ,  Wj  >,  eq.  (3-4) 

maximum  velocity  in  the  boundary  layer,  eq.  (3-40) 

rotation  matrix,  eq.  (3-44) 
mixture  gas  constant 

universal  gas  constant 

A 

species  s  gas  constant  =  R  /  M, 

Reynolds  number 
time 

wedge  thickness.  Demonstration  Case  1 
temperature 

X-,  y-,  and  z-direction  velocity  components 

species  diffusion  velocities  in  the  x-,  y-,  and  z-directions,  eqs.  (3- 
2),  (3-4) 

compressibility  corrected  velocity,  eq.  (3-41) 
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Ux 


U' 


u 

V 

w 

w' 

w 

(x.y.z) 

X. 

y* 


Z' 


z 

Zi 


friction  velocity,  cq.  (3-38) 

vector  of  conserved  variables,  eq.  (3-1) 

transformed  vector  of  conserved  variables,  eqs.  (3-9),  (3-10) 

nonconservative  solution  vector,  eq.  (3-23) 

Cole’s  tabulated  wake  function,  eq.  (3-41) 

vector  of  chemistry  source  terms,  eq.  (3-1) 

transformed  vector  of  chemistry  source  terms,  eqs.  (3-9),  (3-10) 

physical  coordinate  system 

species  s  molar  concentration,  eq.  (3-4) 

law-of-the-wall  coordinate,  eq.  (3-38) 

nontransformed  boundary  condition  vector,  eq.  (3-23) 

transformed  boundary  condition  vector,  eqs.  (3-26),  (3-27) 

curve-fit  coefficients,  eqs.  (4-1)  -  (4-3) 
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Greek 


Oij 

5ij. 

6V 


At 

AU" 

Y 


Ys 


r 

K 

X 

Ub 


Qj,k) 


Jt 


defined  by  eq.  (4-13) 

Kronekar  delta,  eq.  (3-3) 
implicit  delta,  eq.  (3-12) 
nonconservative  flux  vector  delta 
defined  by  eq.  (4-7) 
time  step 

explicit  delta,  eq.  (3-12) 

transformed  coordinate  system 

ratio  of  specific  heats 

species  s  catalytic  efficiency,  eq.  (3-31) 

nontransformed  boundary  condition  vector,  eq.  (3-23) 

von  Karman  constant,  eq.  (3-37) 

thermal  conductivity 

viscosity,  eq.  (3-3) 

bulk  viscosity,  eq.  (3-3) 

collision  cross-section,  eqs.  (4-7),  (4-9),  (4-10) 

constant,  eq.  (3-41) 


84 


Ps 


p 


i 


o 


CO, 

(0 


X 


density  of  species  s 

total  density  =  summation  of  species  densities 
compressibility  correction,  cq.  (3-41) 
shear  stress  tensor,  eq.  (3-3) 
species  s  production/depletion  rate 
vorticity,  eq.  (3-37) 

transformed  boundary  condition  vector,  eqs.  (3-26),  (3-27) 


\ 
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Subscript 

adj 

adjacent  cell 

be 

boundary  cell 

e 

edge 

oo 

freestrcam 

KLEB 

Klebanoff 

max 

maximum 

mix 

mixture 

o 

stagnation  property 

SL 

sea  level 

trans 

transpon 

turb 

turbulent 

VD 

van  Driest 

wall 

wall  value 
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nontransformed  (physical)  variables 


positive  fluxes 
negative  fluxes 
transpose 
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Units 


atm 

ft 

in 

J 

kg 

K 

m 

N 

psia 

R 

sec 


atmosphere 

feet 

inch 

Joule 

kilogram 

Kelvin 

meter 

Newton 

pounds  (force)  per  square  inch 

degrees  Rankine 

second 
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APPENDIX:  THERMOCHEMISTRY  MODEL  COEFFICIENTS 


This  appendix  contains  the  coefficients  for  the  thermochemistry  models  described  in  section  4.0. 
Section  A.  1  contains  the  coefficients  for  the  chemistry  model,  section  A.2  contains  coefficients 
for  the  species  thermodynamics  properties,  and  section  A.3  contains  coefficients  for  the  transport 
properties  equations.  For  details  of  the  format  of  these  tables,  refer  to  the  discussion  in  the 
Software  User’s  manual  (ref.  3). 


A-l 


A.l  Chemistry  Model  Coefficients 


This  section  contains  the  default  chemistry  model  delivered  with  the  HYUDA  code.  The  model 
is  the  modified  Park  and  Menees  model  (modified  by  R.S.  Wong  and  T.R.  Bussing,  see  the 
Interim  Technical  Report,  ref.  2),  with  11  species  (including  electrons)  and  20  reaction  paths.  i 

Reaction  rate  coefficients  are  of  the  Arrhenius  form. 
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AJ  Species  Thermodynamics  Coefficients 


This  section  contains  the  default  species  thermodynamics  database  delivered 
with  the  HYLDA  code.  The  database  curve  fits  are  based  on  the  McBride 
curve  fits. 
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A  J  Transport  Properties  Coefficients 

This  section  contains  the  default  transport  properties  database  delivered 
with  the  HYLDA  code.  The  database  is  based  on  the  Nicolet  transport 
property  model  for  viscosity,  thermal  conductivity,  and  species 
diffusion  coefficients. 

NIOOLET  TRANSPORT  PROPERTY  DATA  (MJ,  LA®DA,  DIFF.  OCEFF.) 
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